      COMMON/CONTRL/NSTEP,NOR(10,3)
      COMMON/VISCOS/CAPPA,KA(5)
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/FLIPS/IFLIPX,IFLIPT,IFLIPR,ICYCLE
      COMMON/UPBOND/PDOWN,UPMACH
      COMMON/PROP/GAMMA,W,G1
      COMMON/ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
      COMMON/CYCLE/ICX,ICT,ICR,IMX,IMT,IMR
      COMMON/PRESS/P(3,18,4),LR,IDLOW,IDHIGH
      COMMON/DIST/TDA(20),RREF(20),UTHETA(20),INOPT
      COMMON/US1/U1(1728,100),UA1(120,100),JAC1(20,100)
      COMMON/UAB/RAD(20),TSET(20),RADIUS(20),DYDX(20),
     1D2YDX(20),S(20),EM(20),WORK(3556)
      REAL JAC1
      REAL NAME(3)
      DATA NAME /'X OP','T OP','R OP'/
C
      IUSE=0
      IMX=1
      IMT=1
      IMR=1
      DO 10 J=1,3
      DO 10 I=1,10
10    NOR(I,J)=0
C
      READ(5,500)(TITLE(I),I=1,80)
500   FORMAT(80A1)
      READ(5,501)NBLADE
501   FORMAT(5I5)
      IF(NBLADE.LE.0)NBLADE=1
      READ(5,501)NX,NTH,NR
      READ(5,501)N1,N2,KUTTA,IEXPLE,IEXPTE
      READ(5,501)IDLOW,IDHIGH,IFCL,IROT
      IF(IFCL.NE.1)IFCL=0
      IF(IROT.NE.0)IROT=1
      IF(IDLOW.LE.0)IDLOW=1
      IF(IDHIGH.LE.0)IDHIGH=NR
      DX=2.0/(N2-N1)
      DR=1.0/(NR-2)
      DTH=1.0/(NTH-2)
      READ(5,502)GAMMA,W,CAPPA,PDOWN
502   FORMAT(8F10.5)
C
      CALL OPEN
C
      READ(5,501)INOPT
      IF(INOPT.NE.1)INOPT=2
      READ(5,501)NUTH
      READ(5,506)(RAD(I),TSET(I),I=1,NUTH)
506   FORMAT(2F10.5)
      DO 11 L=1,NR
      RADIUS(L)=UA1(L,1)
      IF(NUTH.EQ.1)UTHETA(L)=TSET(1)
11    CONTINUE
      IF(NUTH.GT.1)CALL SPLINT(NUTH,NR,RAD,TSET,RADIUS,
     1             UTHETA,DYDX,D2YDX,S,EM,WORK)
      READ(5,501)NTT
      READ(5,506)(RAD(I),TSET(I),I=1,NTT)
      DO 12 L=1,NR
      IF(NTT.EQ.1)TDA(L)=TSET(1)
12    CONTINUE
      IF(NTT.GT.1)CALL SPLINT(NTT,NR,RAD,TSET,RADIUS,
     1            TDA,DYDX,D2YDX,S,EM,WORK)
      READ(5,501)NPT
      READ(5,506)(RAD(I),TSET(I),I=1,NPT)
      DO 13 L=1,NR
      IF(NPT.EQ.1)RREF(L)=TSET(1)
13    CONTINUE
      IF(NPT.GT.1)CALL SPLINT(NPT,NR,RAD,TSET,RADIUS,
     1            RREF,DYDX,D2YDX,S,EM,WORK)
      READ(5,506)UPMACH
C
      DO 15 I=1,5
15    KA(I)=1
      READ(5,503)DT,NSTEP,IBEG
503   FORMAT(F10.5,2I5)
      IF(IBEG.NE.0)IBEG=1
      IF(UPMACH.LE.0)IBEG=0
      DO 20 I=1,10
      READ(5,501)(NOR(I,J),J=1,3)
      IF(NOR(I,1).LE.0)GO TO 25
      IF(NOR(I,3).LE.0)NOR(I,3)=1
20    CONTINUE
25    CONTINUE
C
      G1=GAMMA-1
C
      WRITE(6,315)(TITLE(I),I=1,80),NBLADE
315   FORMAT(' BLADE3D AXIAL COMPRESSOR INVISCID ANALYSIS CODE.',/,
     *80A1/,' THE NUMBER OF BLADES IS ',I5)
      WRITE(6,316)NSTEP,DT
316   FORMAT(/,I5,' CYCLES ARE TO BE RUN WITH DT = ',F10.5)
      IF(IFCL.EQ.0)WRITE(6,340)
      IF(IFCL.NE.0)WRITE(6,341)
340   FORMAT(' IFCL  = 0: NON-CONSERVATIVE FORM OF EQUATIONS')
341   FORMAT(' IFCL NE 0: CONSERVATIVE FORM OF EQUATIONS,'
     1      ,'            NUMERIC GRID METRICS RECOMMENDED')
      IF(IROT.EQ.0)WRITE(6,342)
      IF(IROT.NE.0)WRITE(6,343)
342   FORMAT(' IROT  = 0: UNSTEADY ENERGY EQUATION')
343   FORMAT(' IROT NE 0: CONSTANT ROTHALPY EQUATION')
      WRITE(6,321)
321   FORMAT(' OPERATOR SEQUENCE IS  '
     */' OPERATOR,DT MULTIPLER,REPEAT COUNT')
      DO 100 I=1,10
      IF(NOR(I,1).EQ.0)GO TO 100
      NOR1 = NOR(I,1)
      WRITE(6,317)NAME(NOR1),NOR(I,2),NOR(I,3)
317   FORMAT(4X,A4,10X,I3,10X,I3)
100   CONTINUE
      WRITE(6,318)NX,NTH,NR
318   FORMAT(/,' GRID POINTS IN X DIRECTION IS ',I5/
     *' GRID POINTS IN T DIRECTION IS ',I5/
     *' GRID POINTS IN R DIRECTION IS ',I5)
      WRITE(6,319)N1,N2,KUTTA
319   FORMAT(' THE BLADE LEADING AND TRAILING EDGES ARE',
     *'  AT AXIAL STATIONS ',I5,' AND ',I5/
     *' THE KUTTA POINT IS LOCATED AT AXIAL STATION',I6)
      WRITE(6,322)IEXPLE,IEXPTE
322   FORMAT(' FIRST DAMPER STATION IS AT ',I5/
     *       ' LAST  DAMPER STATION IS AT ',I5)
      WRITE(6,325)IDLOW,IDHIGH
325   FORMAT(' IDLOW  IS ',I5/' IDHIGH IS ',I5)
      WRITE(6,323)GAMMA,W,CAPPA
323   FORMAT(/,' GAMMA  IS ',F10.5/' BLADE SPEED (WRT/AO) IS ',F10.5,
     */      ' ARTIFICAL VISCOSITY PARAMETER IS ',F10.5)
      WRITE(6,324)PDOWN,UPMACH
324   FORMAT(' DOWNSTREAM PRESSURE TO BE SET TO ',F10.5/
     *       ' UPSTREAM MACH NUMBER IS          ',F10.5)
      WRITE(6,334)(L,RADIUS(L),UTHETA(L),TDA(L),RREF(L),L=1,NR)
334   FORMAT(//'     L','  RADIUS  ','  UTHETA  ',
     1         '   TTIN   ','  PTIN OR JPLUS'/
     2       (I5,4F10.5))
      IF(IBEG.EQ.1)WRITE(6,326)
326   FORMAT(/,' A STARTING SOLUTION WILL BE GENERATED'/,)
      IF(DT.EQ.0)WRITE(6,327)
327   FORMAT(/,' DT SET AUTOMATICALLY TO RECOMMENDED DT BELOW',/)
C
      CALL START(IBEG)
C
      CALL MTHREE
C
      STOP
      END
      SUBROUTINE BLADEP(R1,U,I)
C
C       THIS SUBROUTINE CALCULATES THE PRESSURES ON THE BALDE
C       SURFACE
C
      REAL R1(18),U(5,17,18)
C
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /PASS/K,L,R
      COMMON /PRES/P(3,18,4),IDUM1,IDUM2,IDUM3
C
      NTH1=NTH-1
C
      DO 20 KK=1,2
      K=2
      IF(KK.NE.1)K=NTH1
C
      DO 20 L=1,NR
      R=R1(L)
      CALL FINDP(U,P1,V)
      P(I,L,KK)=P1
20    CONTINUE
C
      RETURN
      END
      SUBROUTINE BBR(U,RV,R,IX,RR1,RX1,XX1,XR1,INS,RTAUH,RTAUT)
C
C     BBR CONTROLS THE R-STEP WALL BOUNDARY CONDITIONS
C
      REAL U(5,17,18),RV(4),R(18),INR,INZ
      REAL RR1(18)
      REAL TT(16),TR(15,16),TX(15,16)
      REAL RX1(18),XX1(18),XR1(18)
      REAL CUR(16,2,2),COSN(16,3,2)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
      COMMON /PRESS/P(3,18,4),LR,IDLOW,IDHIGH
      COMMON/DPM/NRLOW,NRHIGH
C
      CALL TRFIL(TT,TR,TX,CUR,COSN,IX)
      INS=IX
10    CONTINUE
C
C     DO INNER WALL
C
      L=NRLOW+1
      LB=L-1
      IF(IFCL.EQ.0)LB=L
      LR=1
      IF(L.NE.2)LR=3
      ITB=0
      INR=RV(1)
      INZ=RV(2)
      XX=XX1(LB)
      XR=XR1(LB)
      RR=RR1(LB)
      RX=RX1(LB)
      CALL BOUNDR(L,U,INR,INZ,R,ITB,IX,RR,RX,TR,TX,XR,XX,RTAUH)
C
C     DO OUTER WALL
C
      L=NRHIGH-1
      LB=L+1
      IF(IFCL.EQ.0)LB=L
      LR=2
      IF(L.NE.(NR-1))LR=4
      ITB=1
      INR=RV(3)
      INZ=RV(4)
      XX=XX1(LB)
      XR=XR1(LB)
      RX=RX1(LB)
      RR=RR1(LB)
      CALL BOUNDR(L,U,INR,INZ,R,ITB,IX,RR,RX,TR,TX,XR,XX,RTAUT)
C
      RETURN
      END
      SUBROUTINE BOUNDR(L2,U,INR,INZ,R,ITB,IX,RR,RX,Q1,Q2,XR,XX,RTAU)
C
C
C              THIS ROUTINE APPLIES THE HARD WALL BOUNDARY
C              CONDITION FOR THE R TIME STEP
C
      REAL U(5,17,20),INR,INZ,IPR,IPZ,R(18)
      REAL Q1(15,16),Q2(15,16)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /PRESS/P(3,18,4),LR,IDLOW,IDHIGH
      COMMON/PASS/K,L,R1
      COMMON/DPM/NRLOW,NRHIGH
      COMMON/PROP/GAMMA,W,G1
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXTPE
      L1=L2-1
      L3=L2+1
      L=L2
C
C              IS IT TOP OR BOTTOM
C
      IF(ITB.GT.0)GO TO 5
C
C              BOTTOM
C
      LNN=L1
      LNP=L3
      LPQ=1
      IF(LR.NE.1)LPQ=19
      ISN=1
      GO TO 6
C
C              TOP
C
5     LNN=L3
      LNP=L1
      LPQ=NR
      IF(LR.NE.2)LPQ=20
      ISN=-1
6     CONTINUE
C
      NT1=NTH-1
      DO 100 KS=2,NT1
      K=KS
      K1=K-1
      K3=K+1
C
C              GET NORMAL AND PARALLEL VELOCITIES
C              AND PARALLEL SURFACE VECTOR,IPR,IPZ
C
      CALL PROJR(U,K,L,INR,INZ,IPR,IPZ,VP,VN)
      VN=-VN
C
C              GET AXIAL AND RADIAL VELOCITIES
C
      CALL UNPROJ(UR,UZ,VP,VN,INR,INZ,IPR,IPZ)
      RHOR=U(1,K,L)
      UTH = U(3,K,L)/RHOR
      R1=R(L)
      UTHDUM=UTH*R1/R(L-ISN)
      UTH=0.5*(UTH+UTHDUM)
      RHO=RHOR/R1
      DPDX=(P(3,K,LR)-P(1,K,LR))/2./DX
      K=K1
      CALL FINDP(U,PK1,V)
      K=K3
      CALL FINDP(U,PK3,V)
      K=KS
      DPDTH=(PK3-PK1)/2./DTH
      RS=R1*ISN
      TEM1=INR*UTH**2/RS
      TEM2=RTAU*VP**2
      DPDN=RHO*(TEM1+TEM2)
      TR=Q1(K1,L-1)
      TX=Q2(K1,L-1)
      DPDX=DPDX*(XR*INR+XX*INZ)
      DPDTH=DPDTH*(TR*INR+TX*INZ)
      T=RR*INR+RX*INZ
      R1=R(LNP)
      CALL FINDP(U,P1,V)
      R1=R(LNN)
      CALL FINDP(U,P3,V)
      R1=R(L)
      CALL FINDP(U,P2,V)
      DAMP=(PK3+PK1-2.*P2) + (P3+P1-2.*P2)
      DAMP=DAMP+ (P(3,K,LR)+P(1,K,LR)-2.*P2)
      DPDR=(DPDN-DPDX-DPDTH-0.4*DAMP)/T
      P1=P2-DPDR*DR*ISN
      URF=AMIN1(P3,1.0)
      P1=P3+URF*(P1-P3)
C
C              SET UP DUMMY POINTS
C
      RS=R1
      R1=R(LNN)
      UTH=UTHDUM
      U(2,K,LNN)=UR
      U(3,K,LNN)=UTH
      U(4,K,LNN)=UZ
      V=UR*UR+UZ*UZ
      V=V+UTH**2
      V1=0.5*V
      RHOM=RHOR/R1
      EM=P1/G1/RHOM
      U(1,K,LNN)=RHOR
      ES=RHOR*(EM+V1)
      U(5,K,LNN)=ES
      DO 10 I=2,4
      U(I,K,LNN)=U(I,K,LNN)*RHOR
10    CONTINUE
C
C
C
100   CONTINUE
      RETURN
      END
      SUBROUTINE BBT(U,INX,R1,TT,TR,TX,CUR,COSN,INS,RR,RX,XX,XR)
C
C     THIS ROUTINE CONTROLS THE THETA BOUNDARY CONDITIONS
C
      REAL U(5,17,18),TT(16),TR(15,16),TX(15,16)
      REAL R1(18)
      REAL RR(18),RX(18),XX(18),XR(18)
      REAL CUR(16,2,2),COSN(16,3,2)
      COMMON/PASSCO/GR,GT,GZ,AR,AT,AZ,BR,BT,BZ
      COMMON /PRES/P(3,18,4),IDUM1,IDUM2,IDUM3
      COMMON/PROP/GAMMA,W,G1
      COMMON/PASS/K,L,R
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/FLIPS/IFLIPX,IFLIPT,IFLIPR,ICYCLE
      NR1=NR-1
      NT1=NTH-1
C
C     IF INX=INS GEOM ALREADY IN CORE
C
      CALL TRFIL(TT,TR,TX,CUR,COSN,INX)
      INS=INX
C
C     IS THIS INSIDE BLADE ROW
C
      IF(INX.GE.N1 .AND. INX.LE.N2)GO TO 100
C
C     UPSTREAM OR DOWNSTREAM OF BLADE ROW, APPLY PERIODICITY
C
      DO 10 L=1,NR
      DO 10 I=1,5
      U(I,1,L)=U(I,NT1,L)
      U(I,NTH,L)=U(I,2,L)
10    CONTINUE
C
C
      RETURN
C
100   CONTINUE
C
      NY=NT1-2
      DO  200 LS=2,NR1
      L=LS
      L1=L-1
      L3=L+1
      LL=L-1
      R=R1(L)
      KK=1
      DO  180 KS=2,NT1,NY
      K=KS
      IF(K.NE.2)KK=2
      K1=K-1
      K3=K+1
      ISN=1
      IF(KS.NE.2)ISN=-1
C
C     FIND DIRECTION COSIN*S OF PARALLEL VECTORS
C     N= BLADE NORMAL VECTOR = > (GR,GT,GZ)
C     A= PARALLEL VECTOR, NEARLY RADIAL , => (AR,AT,AZ)
C     B= PARALLEL VECTOR, RADIAL COMPONENT NEARLY ZERO,=> (BR,BT,BZ)
C
      GR=COSN(LL,1,KK)
      GT=COSN(LL,2,KK)
      GZ=COSN(LL,3,KK)
C
      T= -GR/GT
      AR=1.+ T*T
      AR= SQRT(1./AR)
      AT= T*AR
      AZ=0.
C
      AA= GT-GR/AR*AT
      AA= -GZ/AA
      T= AA*AA
      T1=AT/AR
      BZ=1.+T+T1*T1*T
      BZ= SQRT(1./BZ)
      BT = AA*BZ
      BR=-T1*BT
C
      CALL PROJT(VN,U,GR,GT,GZ)
      CALL PROJT(VPA,U,AR,AT,AZ)
      CALL PROJT(VPB,U,BR,BT,BZ)
      VN = -VN
C
C     SET UP INITIAL GUESS FOR UNPROT
C
      RHOR=U(1,K,L)
      RHO=RHOR/R
      UR=U(2,K,L)/RHOR
      UTH=U(3,K,L)/RHOR-W*R
      UZ=U(4,K,L)/RHOR
      CALL UNPROT(VN,VPA,VPB,UR,UTH,UZ)
      VNT=UR*GR+UTH*GT+UZ*GZ
      KZ =1
      IF(K.NE.2)KZ=NTH
      DPDR=(P(2,L3,KK)-P(2,L1,KK))/2./DR
      DPDX=(P(3,L,KK)-P(1,L,KK))/2./DX
      TEM1=(VPB**2*CUR(LL,1,KK)+VPA**2*CUR(LL,2,KK))
      TEM2=(R*W*W+2.*UTH*W)*GR-2.*UR*W*GT
      DPDN=RHO*(TEM1+TEM2)
      T1=GR*RR(L)+GZ*RX(L)
      DPDR=DPDR*T1
      T3=GR*XR(L)+GZ*XX(L)
      DPDX=DPDX*T3
      T2=GR*TR(K1,LL)+GT/R*TT(LL)+GZ*TX(K1,LL)
      IF(ABS(T2).GT.0.05)GO TO 175
      WRITE(6,600)
600   FORMAT(' TI IS TOO SMALL')
      T2=0.05
175   CONTINUE
      DPDTH=(DPDN-DPDR-DPDX)/T2
      K=KS
      CALL FINDP(U,P2,V)
      P1=P2 -DPDTH*ISN*DTH
      K=KZ
      CALL FINDP(U,P3,VVV)
      TMP=ICYCLE/10.
      URF=AMIN1(P3,1.,TMP)
      IF(P1.LE.0)URF=0.
      UTH=UTH +W*R
      P1=P3+URF*(P1-P3)
      K=KS
      V=UR**2 +UTH**2 +UZ**2
      V1=0.5*V
      EM=P1/G1/RHO
      ES=RHOR*(EM+V1)
      U(1,KZ,L)=RHOR
      U(2,KZ,L)=UR*RHOR
      U(3,KZ,L)=UTH*RHOR
      U(4,KZ,L)=UZ*RHOR
      U(5,KZ,L)=ES
      K=KZ
      CALL FINDP(U,P1,VVV)
180   CONTINUE
C
C
      IF(INX.NE.KUTTA)GO TO 190
      K=2
      CALL FINDP(U,P1,V1)
      K=NT1
      CALL FINDP(U,P2,V2)
      PBAR=0.5*(P1+P2)
      P3=2.0*PBAR-P1
      P4=2.0*PBAR-P2
      WA=W*U(1,1,L)
      V1=(U(2,1,L)**2+(U(3,1,L))**2+U(4,1,L)**2)
      V1=V1/U(1,1,L)
      WA=W*U(1,NTH,L)
      V2=U(2,NTH,L)**2+U(4,NTH,L)**2
      V2=V2+(U(3,NTH,L))**2
      V2=V2/U(1,NTH,L)
      U(5,1,L)=P3/G1*R+V1*0.5
      U(5,NTH,L)=P4/G1*R+V2*0.5
190   CONTINUE
C
200   CONTINUE
      RETURN
      END
      SUBROUTINE XBOUND(IFQ)
C     THIS ROUTINE DOES THE UPSTREAM AND DOWNSTREAM BOUNDARY CONDITIONS
      COMMON /PASS/K,L,R
      COMMON /DIST/TDA(20),RREF(20),UTHETA(20),INOPT
      COMMON /UAB/UA,UB,TU1,TU2
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /IOSAV/IER,IASSOC,IASY
      COMMON /UPBOND/PDOWN,UPMACH
      COMMON /ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
      COMMON /PROP/GAMMA,W,G1
      REAL UA(120),UB(120),JMIN,MACH,JAC1(18),JAC2(18)
      REAL U1(5,17,18),R1(18),RHO0(18),PI(18),JA(18)
      REAL JPLUS,JAC(18)
      REAL TU1(1728),TU2(1728)
      REAL U2(5,17,18),R2(18)
      REAL RR(18),XR(18)
      REAL CX(18),JMIN1,XX(18)
      REAL VTH(18)
      EQUIVALENCE (UA(1),R1(1)),(UA(19),XX(1)),(UA(37),XR(1))
      EQUIVALENCE (UB(1),R2(1)),(UB(73),RR(1))
      EQUIVALENCE (TU1,U1),(TU2,U2)
C
      IR=1
      IW=2
      NTH1 = NTH-1
      NTH2 = NTH-2
      RNTH2= 1.0/FLOAT(NTH2)
      NR1  = NR-1
      NR2  = NR-2
      RNR2 = 1.0/FLOAT(NR2)
      RGAM = 1.0/GAMMA
      RG1  = 1.0/G1
      ALP  = 2.0*RG1
C
      J=NX-1
      CALL IO(IR,UA,IFQ,J,TU1,JAC1)
      CALL IO(IR,UB,IFQ,NX,TU2,JAC2)
      J=NX-2
      IGO=1
      IFW=1
C
5     DO 10 L=1,NR
      PA   =0.0
      TVM  =0.0
      CA   =0.0
      SPEED=0.0
      R=R1(L)
      DO 9 K=2,NTH1
      CALL FINDP(U1,P,V)
      E=(U1(5,K,L)-0.5*V)/U1(1,K,L)
      T=GAMMA*G1*E
      C=SQRT(T)
      PA=PA+P
      SPEED=SPEED+C
      VM=U1(4,K,L)/U1(1,K,L)
      TVM=TVM+VM
      CA=VM+ALP*C*IFW + CA
9     CONTINUE
C     TH AVGS:(PI=PSTATIC, JA=J+(UP) OR J-(DOWN), CX=UX, JAC=A)
      PI(L)=PA*RNTH2
      JA(L)=CA*RNTH2
      CX(L)=TVM*RNTH2
      JAC(L)=SPEED*RNTH2
10    CONTINUE
      CALL PHPROP(U2,R2,0.0)
C
      GO TO (11,30),IGO
C
C     DOWNSTREAM BOUNDARY CONDITIONS
C
C     QUANTITIES FOR TANGENTIAL MOMENTUM: DP/DR=RHO*UTH**2/R
11    CONTINUE
      DO 15 L=2,NR1
      RHO   =0.0
      R     =R1(L)
      V2    =0.0
      VTH(L)=0.0
      DO 14 K=2,NTH1
      V2=U1(3,K,L)**2 +V2
      VTH(L)=VTH(L)+U1(3,K,L)/U1(1,K,L)
      RHO=U1(1,K,L)+RHO
14    CONTINUE
      U2(1,1,L)=RNTH2*RHO/R1(L)
      U2(3,1,L)=NTH2*V2/RHO**2
      VTH(L)   =VTH(L)*RNTH2
15    CONTINUE
C     INTEGRATE TANGENTIAL MOMENTUM
      P=PDOWN
      U2(1,NTH,2)=PDOWN
      DO 20 L=2,NR1
      IF(L.LT.3)GO TO 18
      L1=L-1
      DPDR=0.5*(U2(1,1,L1)*U2(3,1,L1)/R2(L1)+U2(1,1,L)*U2(3,1,L)/R2(L))
      PK=P+(DPDR+PI(L)*XR(L)/DX)*DR/RR(L)
      PK1=1.+DR/DX*XR(L)/RR(L)
      P=PK/PK1
18    U2(1,NTH,L)=P
      PA=P/PI(L)
      RHO=U2(1,1,L)*(PA**RGAM)
      U2(1,2,L)=RHO
      E=P/RHO*RG1
      T=GAMMA*G1*E
      C=SQRT(T)
      UX=JA(L)-ALP*C
      DO 20 K=2,NTH1
      U2(1,K,L)=RHO
      U2(2,K,L)=0.0
      U2(3,K,L)=VTH(L)
      U2(4,K,L)=UX
      U2(5,K,L)=E+0.5*(UX**2+VTH(L)**2)
20    CONTINUE
      CALL CPROP(U2,R2,0.)
      CALL IO(IW,UB,IFQ,NX,TU2,JAC2)
C
C      UPSTREAM BOUNDARY CONDITIONS
C
      CALL IO(IR,UB,IFQ,1,TU2,JAC2)
      IGO=2
      IFW=-1
C     ITERATION COUNTER STORED IN TU2(1728)
      TU2(1728)=TU2(1728)+1
      CALL IO(IR,UA,IFQ,2,TU1,JAC1)
C     FIND SPEED OF SOUND AND UPSTREAM RUNNING CHARACTERISTIC
      GO TO 5
30    CONTINUE
      DO  40 L=2,NR1
      UTH2 = UTHETA(L)**2
      JMIN1=CX(L)-ALP*JAC(L)
      DELX=DT*(JAC(L)-CX(L))
      CAPX=DX/XX(L)
      JMIN=JMIN1+(JA(L)-JMIN1)*(CAPX-DELX)/CAPX
      IF(INOPT.EQ.1)GO TO 190
C
C     JPLUS SPECIFIED BC
C
C     FIND UX AND A FROM INTERSECTION OF CHARACTERISTICS
      JPLUS=ABS(RREF(L))
      UX = 0.5*(JPLUS+JMIN)
      A  = 0.5*(JPLUS-JMIN)/ALP
      V2 = UX**2+UTH2
      RM2= V2/A**2
      T  = TDA(L)/(1.+0.5*G1*RM2)
      GO TO 195
C
C     P0, T0 SPECIFIED BC
C
190   CONTINUE
      AMIN = 0.92
      AMAX = 1.0
      A    = JAC(L)
      ITER = 0
191   ITER = ITER+1
      T=A*A
      JPLUS=JMIN+2.*ALP*A
      UX = 0.5*(JPLUS+JMIN)
      V2 = UX**2+UTH2
      RM2= V2/A**2
      TGUESS=TDA(L)/(1.+0.5*G1*RM2)
      IF(ABS(TGUESS-T).LT.0.0001) GO TO 195
      IF(T.GT.TGUESS)GO TO 193
      AMIN=A
      GO TO 194
193   AMAX=A
194   A=(AMAX+AMIN)*0.5
      IF(ITER.GT.20)GO TO 195
      GO TO 191
C
195   CONTINUE
C     DENSITY FROM ISENTROPIC RELATION
      RHO=T**RG1
C     RESET DENSITY FOR P0 SPECIFIED BC
      IF(INOPT.EQ.1)RHO=RHO*RREF(L)
C     TOTAL ENERGY
      ET = RGAM*RG1*T+0.5*V2
C
      DO 35 K=1,NTH
      U2(1,K,L)=RHO
      U2(2,K,L)=0.0
      U2(3,K,L)=UTHETA(L)
      U2(4,K,L)=UX
      U2(5,K,L)=ET
35    CONTINUE
40    CONTINUE
      CALL CPROP(U2,R2,0.0)
      CALL IO(IW,UB,IFQ,1,TU2,JAC2)
      RETURN
      END


      SUBROUTINE CAFGHK(U,F,G,H,KK,P,V)
      REAL U(5,17,18)
      REAL F(5),G(5),H(5),KK(5)
      COMMON /PROP/GAMMA,WSAV,G1
      COMMON /ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
      COMMON /PASS/K,L,R
      Q =1./R
      U1=U(1,K,L)
      U2=U(2,K,L)
      U3=U(3,K,L)
      U4=U(4,K,L)
      U5=U(5,K,L)
      RS=1./U1
      V =(U2*U2+U3*U3+U4*U4)*RS
      V2=0.5*V
      UTH=U3*RS
      IF(IROT.NE.0) GO TO 10
      P=(U5-V2)*Q*G1
      GO TO 15
10    CONTINUE
      VR=V2*RS
      E =(1./G1+R*WSAV*UTH-VR)/GAMMA
      P =U1*Q*E*G1
      U5=U1*(E+VR)
15    CONTINUE
      PR = P*R
      UR=U2*RS
      F(1) = U2
      F(2) = U2*UR + PR
      F(3) = U3*UR
      F(4) = U4*UR
      F(5) = UR*(U5+ PR)
      UZ = U4*RS
      H(1) = U4
      H(2) = U2*UZ
      H(3) = U3*UZ
      H(4) = U4*UZ + PR
      H(5) = UZ*(U5 + PR)
      UTH = UTH*Q
      KK(1)=0.0
      KK(2) = U3*UTH + P
      KK(3) = -U2*UTH
      KK(4)=0.0
      KK(5)=0.0
      UTH=UTH-WSAV
      G(1) = U1*UTH
      G(2) = U2*UTH
      G(3) = U3*UTH + P
      G(4) = U4*UTH
      G(5) = (U5 + PR)*UTH+WSAV*PR
      RETURN
      END
      SUBROUTINE CLOSE
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /US1/U1(1728,100),UA1(120,100),JAC1(20,100)
      COMMON /FLIPS/IFLIPX,IFLIPT,IFLIPR,ICYCLE
      COMMON /PROP/GAMMA,W,G1
      COMMON /UPBOND/PDOWN,UPMACH
      DATA AZERO/0.0/
      IUNIT = 7
      REWIND IUNIT
      WRITE(6,631)IUNIT
631   FORMAT(///,'  FINAL SOLUTION WRITTEN TO UNIT',I5)
      WRITE(IUNIT)GAMMA,W,G1,UPMACH,AZERO,PDOWN
      DO 100 J=1,NX
      WRITE(IUNIT)(U1(I,J),I=1,1728)
100   CONTINUE
      WRITE(IUNIT)IFLIPX,IFLIPT,IFLIPR,ICYCLE
      RETURN
      END
      SUBROUTINE CPROP(U,R,W)
C     CONVERT PHYSICAL PROPERTIES TO CONSERVATION FORM
C
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      REAL U(5,17,18),R(18)
C
      DO 10 L=1,NR
      DO 10 K=1,NTH
      RHOR=U(1,K,L)*R(L)
      U(1,K,L)=RHOR
      U(2,K,L)=U(2,K,L)*RHOR
      U(3,K,L)=(U(3,K,L)+W*R(L))*RHOR
      U(4,K,L)=U(4,K,L)*RHOR
      U(5,K,L)=U(5,K,L)*RHOR
10    CONTINUE
      RETURN
      END
      SUBROUTINE FINDP(U,P,V)
C       THIS ROUTINE CALCULATES PRESSURE FROM CONSERVATION VARIABLES
      REAL U(5,17,18)
      COMMON /PROP/GAMMA,WSAV,G1
      COMMON/ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      COMMON /PASS/K,L,R
C       CALCULATE  SQUARED VELOCITY V=(UR**Z + UTH**2 + UZ**2)*RHO*R
      V = (U(2,K,L)**2+U(3,K,L)**2+U(4,K,L)**2)/U(1,K,L)
      IF(IROT.NE.0) GO TO 10
      P = (U(5,K,L)-0.5*V)/R*G1
      RETURN
C
C       CALCULATE P USING CONSTANT RHOTHAPLY
C       FAR UPSTREAM ENTHALPY IS 1/(GAMMA-1)
C
10    CONTINUE
      RHOR=U(1,K,L)
      RHIN = 1.0/RHOR
      UTH=U(3,K,L)*RHIN
      U2=V*0.5*RHIN
      E=(1./G1+R*WSAV*UTH-U2)/GAMMA
      P=RHOR/R*E*G1
      U(5,K,L)=RHOR*(E+U2)
      RETURN
      END
      SUBROUTINE IO(IRW,U,IUORUT,J,V,JAC)
C
C     THIS ROUTINE DOES ALL DISK I/O FOR PROBLEM MATRIX
C     IRW=1 => READ   IRW= 2 => WRITE
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /US1/U1(1728,100),UA1(120,100),JAC1(20,100)
      REAL JAC1
      REAL V(1728),JAC(50)
      REAL U(120)
      IF(IRW.EQ.2) GO TO 100
      DO  10 I=1,1728
10    V(I)=U1(I,J)
      DO  20 I=1,120
20    U(I)=UA1(I,J)
      DO  30 I=1,NR
30    JAC(I)=JAC1(I,J)
C
      RETURN
100   CONTINUE
      DO  200 I=1,1728
200   U1(I,J)=V(I)
      RETURN
      END
      SUBROUTINE BLOCKR
C
C              THIS ROUTINE DOES THE DAMPER CALCULATIONS
C              AND CONTROLS THE R-TIME STEP INTEG.
C
C       THIS VERSION FLIPS THE FORWARD-BACKWARD SEQUENCE
C
C
      REAL UA(120),U1(5,17,20),R1(18),RX1(18),RV1(4)
      REAL RVPASS(4),RVD(4)
      REAL XX1(18),XR1(18)
      REAL UB(120),TUSAVE(1728)
      REAL RK(5)
      REAL TU1(1728),RR1(18)
      REAL UT(5,17,18),JAC1(18),JAC2(18)
      COMMON /IOSAV/IER,IASSOC,IASY
      COMMON /FLIPS/IFLIPX,IFLIPT,IFLIPR,ICYCLE
      COMMON /PRESS/P(3,18,4),LR,IDLOW,IDHIGH
      COMMON/PASS/K,L,R
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/DPM/NRLOW,NRHIGH
      COMMON/FLIP/IFLIP
C
C
C
      EQUIVALENCE (UA(1),R1(1)),(TU1,U1),(UA(55),RX1(1))
      EQUIVALENCE (UA(95),RTAUH),(UA(96),RTAUT)
      EQUIVALENCE (UA(19),XX1(1)),(UA(37),XR1(1))
      EQUIVALENCE (UA(73),RR1),(UA(91),RV1(1))
      EQUIVALENCE(UA(97),IDAMP),(UA(100),RVD)
      EQUIVALENCE (UA(104),RTLOW),(UA(105),RTHIGH)
C
C
C
      IR=1
      IW=2
      RK(1)=0.
      RK(2)=1.
      RK(3)=0.
      RK(4)=0.
      RK(5)=0.
      IER=0
      INS=0
      NR1=NR-1
      NY=NTH-1
      IPX=NX-1
      ITOP=NR+2
      ITOP1=NR+1
C
C       FIND WALL PRESSURES AND STORE
C       RETRIEVE DERIVATIVE SWITCHING INFORMATION
C
      I=1
      CALL IO(IR,UA,0,I,TU1,JAC1)
C
      IFLIP=IFLIPR
      IFLIPR=-IFLIPR
      CALL IO(IW,UA,0,I,TU1,JAC1)
C
      CALL WALLP(UA,TU1,2)
      I=2
      CALL IO(IR,UB,0,I,TUSAVE,JAC2)
      CALL WALLP(UB,TUSAVE,3)
C
C
      DO 500 INX=2,IPX
      DO 20 I=1,120
      UA(I)=UB(I)
20    CONTINUE
      DO 30 I=1,18
      JAC1(I)=JAC2(I)
30    CONTINUE
      DO 40 I=1,1728
40    TU1(I)=TUSAVE(I)
      DO  50 K=1,NTH
      DO 50 J=1,4
      P(1,K,J)=P(2,K,J)
50    P(2,K,J)=P(3,K,J)
      IXR=INX+1
      CALL IO(IR,UB,0,IXR,TUSAVE,JAC2)
      CALL WALLP(UB,TUSAVE,3)
C
C       DO FORWARD AND BACKWARD MACCORMACK STEP IN ORDER
C       SPECIFIED BY IFLIP
C
      IFW=0
      NRLOW=1
      NRHIGH=NR
      IF(IDAMP.EQ.1)GO TO 300
      CALL BBR(U1,RV1,R1,INX,RR1,RX1,XX1,XR1,INS,RTAUH,RTAUT)
      CALL STEPR(U1,R1,RR1,RX1,IFW,INX,UT,RK)
C
      GO TO 490
C
C
300   CONTINUE
      NRLOW=1
      NRHIGH=IDHIGH
      RVPASS(1)=RV1(1)
      RVPASS(2)=RV1(2)
      RVPASS(3)=RVD(3)
      RVPASS(4)=RVD(4)
      CALL BBR(U1,RVPASS,R1,INX,RR1,RX1,XX1,XR1,INS,RTAUH,RTLOW)
      NRLOW=IDLOW
      NRHIGH=NR
      RVPASS(1)=RVD(1)
      RVPASS(2)=RVD(2)
      RVPASS(3)=RV1(3)
      RVPASS(4)=RV1(4)
      CALL BBR(U1,RVPASS,R1,INX,RR1,RX1,XX1,XR1,INS,RTHIGH,RTAUT)
      DO 310 I=1,5
      DO 310 K=1,NTH
      T=U1(I,K,IDHIGH)
      U1(I,K,IDHIGH)=U1(I,K,ITOP)
      U1(I,K,ITOP)=T
310   CONTINUE
      NRLOW=1
      NRHIGH=IDHIGH
      CALL STEPR(U1,R1,RR1,RX1,IFW,INX,UT,RK)
      DO 320 I=1,5
      DO 320 K=1,NTH
      T=U1(I,K,ITOP)
      U1(I,K,ITOP)=U1(I,K,NRHIGH)
      U1(I,K,NRHIGH)=T
      T=U1(I,K,ITOP1)
      U1(I,K,ITOP1)=U1(I,K,IDLOW)
      U1(I,K,IDLOW)=T
320   CONTINUE
      NRLOW=IDLOW
      NRHIGH=NR
      CALL STEPR(U1,R1,RR1,RX1,IFW,INX,UT,RK)
      DO 330 I=1,5
      DO 330 K=1,NTH
      T=U1(I,K,ITOP1)
      U1(I,K,ITOP1)=U1(I,K,IDLOW)
      U1(I,K,IDLOW)=T
330   CONTINUE
490   CONTINUE
      CALL IO(IW,UA,0,INX,TU1,JAC1)
C
C
C
500   CONTINUE
      RETURN
      END
      SUBROUTINE BLOCKT
C
C       THIS ROUTINE CONTROLS THE THETA TIME STEP
C
C       THIS VERSION DOES FORWARD-BACKWARD FLIPPING
C       STEPT IS DONE ON ONE SHOT AND REPEATED CYCLES ARE ALLOWED
C       BLADE PRESSURE CALCULATION IS IMPROVED. SMALLER ARRAY RE-
C       QUIRED.
C
      REAL UA(120),U1(5,17,18),R1(18),TT1(16),TR1(15,16),TX1(15,16)
      REAL UT(5,17,18),TK(5),CUR(16,2,2),COSN(16,3,2)
      REAL TU1(1728),TUSAVE(1728),UB(120)
      REAL XX1(18),XR1(18),RR1(18),RX1(18)
      REAL JAC1(18),JAC2(18)
C
      COMMON/IOSAV/IER,IASSOC,IASY
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON /PASS/K,L,R
      COMMON /PRES/P(3,18,4),IDUM1,IDUM2,IDUM3
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /FLIPS/IFLIPX,IFLIPT,IFLIPR,ICYCLE
      COMMON/FLIP/IFLIP
      COMMON /CYCLE/NCYX,NCYT,NCYR,IMX,IMT,IMR
C
      EQUIVALENCE (UA,R1),(TU1,U1)
      EQUIVALENCE (UA(19),XX1)
      EQUIVALENCE(UA(55),RX1),(UA(37),XR1),(UA(73),RR1)
C
      IR=1
      IW=2
      TK(1)=0.
      TK(2)=0.
      TK(3)=1.
      TK(4)=0.
      TK(5)=0.
      IER=0
      IPX=NX-1
      INXS=0
      NR1=NR-1
      NTH1=NTH-1
      NY=NTH1-2
      IFW=0
      N11=N1-1
      N22=N2+1
      N12=N1-2
      ISP=0
C
C        RETRIEVE DERIVATIVE SWITCHING INFORMATION
C
      IFLIP=IFLIPT
      IFLIPT=-IFLIPT
C
C
      CALL IO(IR,UB,0,2,TUSAVE,JAC2)
C
C       START CYCLE.  DO FORWARD AND BACKWARD SETPS IN ONE SHOT
C
      DO 500 INX=2,IPX
      INX1=INX+1
C
      DO 20 I=1,1728
      TU1(I)=TUSAVE(I)
20    CONTINUE
      DO 30 I=1,120
      UA(I)=UB(I)
30    CONTINUE
      DO 40 I=1,18
      JAC1(I)=JAC2(I)
40    CONTINUE
C
      CALL IO(IR,UB,0,INX1,TUSAVE,JAC2)
C
C
C       CALCULATE BLADE SURFACE PRESSURES WHEN INSIDE BLADE ROW.
C
      IF(INX.LT.N12 .OR. INX.GT.N2)GO TO 60
      DO 50 L=1,18
      DO 50 KK=1,2
      P(1,L,KK)=P(2,L,KK)
      P(2,L,KK)=P(3,L,KK)
50    CONTINUE
      CALL BLADEP(UB,TUSAVE,3)
      ISP=1
60    CONTINUE
C
      DO 300 NC=1,NCYT
C
C       UPDATE SURFACE PRESSURES WHEN REPEATING CYCLES
C
      IF(NC.NE.1.AND.ISP.EQ.1)CALL BLADEP(UA,TU1,2)
C
      CALL BBT(U1,INX,R1,TT1,TR1,TX1,CUR,COSN,INXS,
     *       RR1,RX1,XX1,XR1)
C
      DO 250 L=1,NR
      DO 250 K=1,NTH
      DO 250 I=1,5
      UT(I,K,L)=U1(I,K,L)
250   CONTINUE
C
      CALL STEPT(U1,UT,INX,IFW,TT1,TR1,TX1,R1,TK)
C
300   CONTINUE
C
      CALL IO(IW,UA,0,INX,TU1,JAC1)
      ISP=0
C
C
500   CONTINUE
      RETURN
      END
      SUBROUTINE BLOCKX
C
C       THIS ROUTINE CONTROLS THE X-TIME STEP INTEGRATION
C
C       DERIVATIVE SWITCHING IS DONE
C
C       THIS IS THE VERSION THAT WORKS WITH THE NEW VERSION OF
C       STEPX WHICH DOES NOT NEED ANY PROVISIONAL FILES. ARTI-
C       FICIAL VISCOSITIES ARE TREATED AS USUAL.
C
      REAL UA(120),UB(120),UC(120)
      REAL JAC1(18),JAC2(18),JAC3(18)
      REAL TU1(1728),TU2(1728),TU3(1728)
C
      COMMON/UAB/UA,UB,TU1,TU2
      COMMON /IOSAV/IER,IASSOC,IASY
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/FLIP/IFLIP
      COMMON /FLIPS/IFLIPX,IFLIPT,IFLIPR
C
      REAL U1(5,17,18),R1(18)
      REAL U2(5,17,18),R2(18),XX(18),XR(18)
      REAL U3(5,17,18),R3(18)
      REAL UP(5,17,18)
C
      EQUIVALENCE (UA(1),R1(1)),(TU1,U1)
     *       ,(UB(1),R2(1)),(UB(19),XX(1)),(TU2,U2),(UB(37),XR(1)),
     *       (UC(1),R3(1)),(TU3,U3)
      REAL XX1(18),XR1(18)
      EQUIVALENCE (UA(19),XX1(1)),(UA(37),XR1(1))
      REAL XX3(18),XR3(18)
      EQUIVALENCE (UC(19),XX3(1)),(UC(37),XR3(1))
C
      IR=1
      IW=2
      IER=0
C       SET UPSTREAM AND DOWNSTREAM BOUNDARY CONDITIONS
C
      CALL XBOUND
C
      IF(DT.EQ.0)GO TO 1000
C
C       RETRIEVE DERIVATIVE SWITCHING INFORMATION
C
      IFLIP=IFLIPX
      IFLIPX=-IFLIPX
C
      IF(IFLIP.EQ.-1)GO TO 25
C
C       FOR IFLIP=+1 MOVE FROM J = 1  TO  J = NX.
C
      INX=2
      IFW=0
      CALL IO(IR,UA,IFW,1,TU1,JAC1)
      CALL IO(IR,UB,IFW,2,TU2,JAC2)
      CALL IO(IR,UC,IFW,3,TU3,JAC3)
C
C       COPY BOUNDARY CONDS. INTO PROVISIONAL VALUES.
C
      DO 5 L=1,NR
      DO 5 K=1,NTH
      DO 5 I=1,5
      UP(I,K,L)=U1(I,K,L)
5     CONTINUE
C
15    CALL STEPX(U1,U2,U3,UP,R1,R2,R3,XR,XX,IFW,INX)
      IPX=INX+1
      IMX=INX-1
      CALL IO(IW,UB,IFW,IMX,TU1,JAC2)
      IF(IPX.GT.NX)GO TO 20
      CALL IO(IR,UA,IFW,IPX,TU1,JAC1)
      CALL STEPX(U2,U3,U1,UP,R2,R3,R1,XR3,XX3,IFW,INX)
      IPX=INX+1
      IMX=INX-1
      CALL IO(IW,UC,IFW,IMX,TU2,JAC2)
      IF(IPX.GT.NX)GO TO 20
      CALL IO(IR,UB,IFW,IPX,TU2,JAC2)
      CALL STEPX(U3,U1,U2,UP,R3,R1,R2,XR1,XX1,IFW,INX)
      IPX=INX+1
      IMX=INX-1
      CALL IO(IW,UA,IFW,IMX,TU3,JAC3)
      IF(IPX.GT.NX)GO TO 20
      CALL IO(IR,UC,IFW,IPX,TU3,JAC3)
      GO TO 15
C
C       FOR IFLIP=-1 MOVE FROM J = NX  TO  J = 1.
C
25    CONTINUE
      INX=NX-1
      IFW=0
      CALL IO(IR,UA,IFW,NX,TU1,JAC1)
      CALL IO(IR,UB,IFW,NX-1,TU2,JAC2)
      CALL IO(IR,UC,IFW,NX-2,TU3,JAC3)
C
C       WRITE BOUNDARY CONDS. INTO PROVISIONAL VALUE.
C
      DO 6 L=1,NR
      DO 6 K=1,NTH
      DO 6 I=1,5
      UP(I,K,L)=U1(I,K,L)
6     CONTINUE
C
16    CALL STEPX(U1,U2,U3,UP,R1,R2,R3,XR,XX,IFW,INX)
      IPX=INX+1
      IMX=INX-1
      CALL IO(IW,UB,IFW,IPX,TU1,JAC1)
      IF(IMX.LT.1)GO TO 20
      CALL IO(IR,UA,IFW,IMX,TU1,JAC1)
      CALL STEPX(U2,U3,U1,UP,R2,R3,R1,XR3,XX3,IFW,INX)
      IPX=INX+1
      IMX=INX-1
      CALL IO(IW,UC,IFW,IPX,TU2,JAC2)
      IF(IMX.LT.1)GO TO 20
      CALL IO(IR,UB,IFW,IMX,TU2,JAC2)
      CALL STEPX(U3,U1,U2,UP,R3,R1,R2,XR1,XX1,IFW,INX)
      IPX=INX+1
      IMX=INX-1
      CALL IO(IW,UA,IFW,IPX,TU3,JAC3)
      IF(IMX.LT.1)GO TO 20
      CALL IO(IR,UC,IFW,IMX,TU3,JAC3)
      GO TO 16
C
20    CONTINUE
1000  CONTINUE
      RETURN
      END
      SUBROUTINE MD(MDOT,U,R,RR1,RX1,XX1,XR1,JAC,JKK)
      REAL U(5,17,18),MDOT,R(18),DP(18)
      REAL RHO(18),MTH(17)
      REAL JAC(18)
      REAL RR1(18),RX1(18),XX1(18),XR1(18)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      NR1=NR-1
      NTH1=NTH-1
      NTH2=NTH-2
      MDOT=0.
      NR2=NR-2
      IF(JKK.EQ.0)GO TO 91
      IF(JKK.NE.2)GO TO 5
      DO 4 I=1,5
      DO 3 K=1,NTH
      U(I,K,1)=U(I,K,2)
      U(I,K,NR)=U(I,K,NR1)
3     CONTINUE
      DO 4 L=1,NR
      U(I,1,L)=U(I,NTH1,L)
      U(I,NTH,L)=U(I,2,L)
4     CONTINUE
5     CONTINUE
91    CONTINUE
      DO 20 L=1,NR
      DP1=0.
      SINPHI=XR1(L)
      COSPHI=XX1(L)
      IF(L.GT.1.AND.L.LT.NR)GO TO 1091
      L1=2
      IF(L.EQ.NR)L1=NR1
      SINPHI=0.5*(SINPHI+XR1(L1))
      COSPHI=0.5*(COSPHI+XX1(L1))
1091  CONTINUE
      RHO1=0.
      DO 10 K=1,NTH
      KK=K
      RHO1=U(1,K,L)/R(L)
      UX=U(4,K,L)/U(1,K,L)
      UR=U(2,K,L)/U(1,K,L)
      IF(L.GT.1.AND.L.LT.NR)GO TO 1092
      RHO1=U(1,K,L1)/R(L1)
      IF(K.EQ.1)UX=U(4,2,L)/U(1,2,L)
      IF(K.EQ.1)UR=U(2,2,L)/U(1,2,L)
      IF(K.EQ.NTH)UX=U(4,NTH1,L)/U(1,NTH1,L)
      IF(K.EQ.NTH)UR=U(4,NTH1,L)/U(1,NTH1,L)
      UX=0.5*(UX+U(4,K,L)/U(1,K,L))
      UR=0.5*(UR+U(2,K,L)/U(1,K,L))
1092  T=RHO1*R(L)*(UX*COSPHI+UR*SINPHI)
      MTH(K)=T
      RHO1=RHO1+U(1,K,L)
10    CONTINUE
      MTH(1)=MTH(1)*0.5
      MTH(NTH)=MTH(NTH)*0.5
      DO 11 K=2,NTH
      KK=K-1
      DP1=DP1+(MTH(K)+MTH(KK))*0.5
11    CONTINUE
      D=JAC(L)
      IF(L.EQ.1)D=0.5*(D+JAC(2))
      IF(L.EQ.NR)D=0.5*(D+JAC(NR1))
      DP(L)=DP1*DTH
      DP(L)=DP(L)/D
20    CONTINUE
      DP(1)=DP(1)*0.5
      DP(NR)=DP(NR)*0.5
      DO 30 L=2,NR
      L1=L-1
      T=1.0
      MDOT=MDOT+0.5*(DP(L1)+DP(L))*T
30    CONTINUE
C      MDOT CALCULATED AT STP: RHO0=0.07647 LB/FT**3, A0=1116.43 FT/SEC
      MDOT=MDOT*FLOAT(NBLADE)*DR*0.07647*1116.43
      RETURN
      END
      SUBROUTINE MFLOW
C     USES PTASS TO FIND MASS FLOWS THROUGHOUT FIELD
C     AND CALCULATES MEAN AND STANDARD DEVIATION.
      COMMON/BLADE/N1,N2,KUTTA,IEL,IET
      COMMON/GRID/ NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/MASS/ MDSAV(100)
      REAL MDAV,MDSQ,MDUP,MDOA,MDSAV
C     FIND MASS FLOWS
      CALL PTASS(0,-3,1,NX,2,2,0,0)
C     FIND MEAN AND STANDARD DEVIATION
      N1M1=N1-1
      MDAV=0.0
      MDSQ=0.0
C     FIRST UPSTREAM
      DO 10 J=1,N1M1
      MDAV=MDAV+MDSAV(J)
10    MDSQ=MDSQ+MDSAV(J)**2
      MDUP=MDAV/N1M1
      SDUP=SQRT((MDSQ-N1M1*MDUP**2)/(N1-2))
C     THEN OVERALL
      DO 11 J=N1,NX
      MDAV=MDAV+MDSAV(J)
11    MDSQ=MDSQ+MDSAV(J)**2
      MDOA=MDAV/NX
      SDOA=SQRT((MDSQ-  NX*MDOA**2)/(NX-1))
C     OUTPUT
      WRITE(6,601)
      WRITE(6,602)(J,MDSAV(J),J=1,NX)
      WRITE(6,603)MDUP,SDUP
      WRITE(6,604)MDOA,SDOA
601   FORMAT(//,T2,'WHEEL MASS FLOW (LBM/SEC AT STP) AT AXIAL STATION J'
     1,/,T5,5('J      MDOT',4X))
602   FORMAT(5(I5,F10.5))
603   FORMAT(/,' UPSTREAM :     MEAN= ',F10.5,'; STANDARD DEVIATION= ',
     1F10.5)
604   FORMAT(/,' OVERALL  :     MEAN= ',F10.5,'; STANDARD DEVIATION= ',
     1F10.5)
      RETURN
      END
      SUBROUTINE MTHREE
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/RESIDS/RESID(18)
      COMMON /CYCLE/ICX,ICT,ICR,IMX,IMT,IMR
      INTEGER CODE(50),JST(50),JEND(50),LST(50),LEND(50)
      COMMON /BLADE/N1,N2,KUTTA,IEL,IET
      COMMON /CONTRL/NSTEP,NOR(10,3)
      COMMON/FLIPS/IFLIPX,IFLIPT,IFLIPR,ICYCLE
      COMMON /MASS/MDSAV(100)
      REAL MDSAV
      REAL XPLOT(100),RESSAV(100)
      COMMON/PROUT/PRPS(100),PRSS(100)
      INTEGER IAXNUM(17),MTITLE(15)
      INTEGER RESTIT(18),PRESTL(14),CYCLE(10)
      DATA IAXNUM/'A','X','I','A','L',' ','G','R','I','D',
     1            ' ','P','L','A','N','E',0/
      DATA MTITLE/'M','A','S','S',' ','F','L','O','W',
     1            ' ','R','A','T','E',0/
      DATA RESTIT/'T','I','P',' ','T','.','E','.',' ',
     1            'P','R','E','S','S','U','R','E',0/
      DATA PRESTL/'S','T','A','I','C',' ','P','R','E','S',
     1            'U','R','E',0/
      DATA CYCLE /'I','T','E','R','A','T','I','O','N',0/
      DTSAV=DT
      IFLIPX=1
      IFLIPT=1
      IFLIPR=1
      DO 5 J=1,100
5     XPLOT(J)=J
      READ(5,500)NP
      READ(5,500)(CODE(I),JST(I),JEND(I),LST(I),LEND(I),I=1,NP)
      DO 10 I=1,NP
10    CALL PTASS(0,CODE(I),JST(I),JEND(I),LST(I),LEND(I),0,0)
C.... WRITE MASS FLOWS FOR J = 1,NX.
      CALL MFLOW
      READ(5,500)NP,NTIMES
      IF(NP.GT.0)READ(5,500)(CODE(I),JST(I),JEND(I),
     1           LST(I),LEND(I),I=1,NP)
      IF(NTIMES.EQ.0)NP=0
      IF(DT.LE.0)GO TO 210
      ICYCLE = 0
      WRITE(6,610)NSTEP
      WRITE(6,620)(L,L=1,NR)
C
C
      ISC=0
      ISCMAX=100
      IF(ISCMAX.GT.NTIMES.AND.NTIMES.NE.0)ISCMAX=NTIMES
C
C
      IF(NSTEP.LE.0)GO TO 210
      DO 100 N=1,NSTEP
      ISC=ISC+1
      DO 90 I=1,10
      IXRT=NOR(I,1)
      IF(IXRT.EQ.0) GO TO 90
      IR=NOR(I,2)
      IF(IXRT.NE.1) GO TO 20
      IMX=IR
      ICX=NOR(I,3)
      TAUX=DTSAV/DX*IMX
      DT=DTSAV*IMX
      DO 11 II=1,ICX
      ICYCLE=ICYCLE+1
      CALL BLOCKX
11    CONTINUE
C
20    CONTINUE
      IF(IXRT.NE.2)GO TO 30
      IMT=IR
      ICT=NOR(I,3)
      TAUT=DTSAV/DTH*IMT
      DT=DTSAV*IMT
      CALL BLOCKT
C
C
30    CONTINUE
      IF(IXRT.NE.3)GO TO 80
      IMR=IR
      ICR=NOR(I,3)
      TAUR=DTSAV/DR*IMR
      DT=DTSAV*IMR
      DO 13 II=1,ICR
      CALL BLOCKR
13    CONTINUE
80    CONTINUE
90    CONTINUE
      WRITE(6,630)N,(RESID(L),L=1,NR)
      RESSAV(ISC)=RESID(NR-1)
      IF(ISC.LT.ISCMAX)GO TO 93
      CALL PPSETU(ISC,XPLOT,RESSAV,'R',CYCLE,RESTIT,1)
      CALL PTASS(0,-3,1,NX,2,2,0,0)
      CALL PPSETU(NX,XPLOT,PRPS,'P',IAXNUM,PRESTL,-1)
      CALL PPSETU(NX,XPLOT,PRSS,'S',IAXNUM,PRESTL,0)
      CALL PPSETU(NX,XPLOT,MDSAV,'M',IAXNUM,MTITLE,1)
      ISC=0
      IF(N.NE.NSTEP)WRITE(6,620)(L,L=1,NR)
93    CONTINUE
      IF(NTIMES.EQ.0) GO TO 100
      IF(NP.EQ.0)GO TO 100
      IT=N/NTIMES*NTIMES
      IF(N.NE.IT) GO TO 100
      CALL CLOSE
      DO 95 I=1,NP
95    CALL PTASS(0,CODE(I),JST(I),JEND(I),LST(I),LEND(I),0,0)
      WRITE(6,620)(L,L=1,NR)
100   CONTINUE
      WRITE(6,600)NSTEP
      READ(5,500)NP
      READ(5,500)(CODE(I),JST(I),JEND(I),LST(I),LEND(I),I=1,NP)
      DO 200 I=1,NP
200   CALL PTASS(0,CODE(I),JST(I),JEND(I),LST(I),LEND(I),0,0)
210   CONTINUE
C.... WRITE MASS FLOWS FOR J = 1,NX.
      CALL MFLOW
      CALL PPSETU(NX,XPLOT,PRPS,'P',IAXNUM,PRESTL,-1)
      CALL PPSETU(NX,XPLOT,PRSS,'S',IAXNUM,PRESTL,0)
      CALL PPSETU(NX,XPLOT,MDSAV,'M',IAXNUM,MTITLE,1)
      CALL TPM
      CALL CLOSE
500   FORMAT(5I5)
600   FORMAT(///,I5,' OPERATOR CYCLES HAVE BEEN RUN.')
610   FORMAT(///,I5,' OPERATOR CYCLES ARE TO BE RUN.')
620   FORMAT(1H1,T2,'TRAILING EDGE STATIC PRESSURES FROM HUB TO TIP.',
     1        /,1X,T2,'CYC L = ',18(I2,5X))
630   FORMAT(1X,T1,I4,T6,18(F6.4,1X))
      RETURN
      END
      SUBROUTINE OPEN
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /US1/U1(1728,100),UA1(120,100),JAC1(20,100)
      COMMON/TGSAV/TTSAV(18,18,100),TRSAV(18,18,100),TXSAV(18,18,100),
     *CURSAV(18,2,2,100),COSNSV(18,3,2,100)
      REAL TU(1728),UA(150),JAC(50)
      REAL JAC1
      INTEGER XXS,XXE,XRS,XRE,RXS,RXE,RRS,RRE
      DATA XXS,XRS,RXS,RRS/19,37,55,73/
      XXE=NR+18
      XRE=NR+36
      RXE=NR+54
      RRE=NR+72
      DO 100 J=1,NX
      READ(2)((TTSAV(K,L,J),K=1,NTH),L=1,NR),
     *((TRSAV(K,L,J),TXSAV(K,L,J),K=1,NTH),L=1,NR),
     *(((CURSAV(L,KX,KY,J),KX=1,2),(COSNSV(L,KX,KY,J),KX=1,3),
     *KY=1,2),L=1,NR)
      READ(3)(UA1(I,J),I=1,NR),(UA1(I,J),I=XXS,XXE),
     *(UA1(I,J),I=XRS,XRE),(UA1(I,J),I=RXS,RXE),
     *(UA1(I,J),I=RRS,RRE)
      READ(3)(UA1(I,J),I=91,105)
      READ(3)(JAC1(L,J),L=1,NR)
100   CONTINUE
      RETURN
      END
      SUBROUTINE PHPROP(U,R,W)
C
C     CONVERT CONSERVATION VARIABLES TO PHYSICAL VARIABLES
C     E.G. R*RHO*E  =>  E
C
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      REAL U(5,17,18),R(18)
C
      DO 10 L=1,NR
      RR=1.0/R(L)
      DO 10 K=1,NTH
      RHORI=1.0/U(1,K,L)
      U(1,K,L)=U(1,K,L)*RR
      U(2,K,L)=U(2,K,L)*RHORI
      U(3,K,L)=U(3,K,L)*RHORI-W*R(L)
      U(4,K,L)=U(4,K,L)*RHORI
      U(5,K,L)=U(5,K,L)*RHORI
10    CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(U,IMOD,MSG,J,LS,LSS,DMD,R1)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/PROP/GAMMA,W,G1
      REAL U(5,17,18),MSG(5,5)
      REAL PTS(18),TTS(18)
      REAL R1(18)
      REAL R(18)
      REAL AV(5)
      TODEG=57.29578
      G2=G1/GAMMA
      NR1=NR-1
      NTH1=NTH-1
      NR2=NR-2
      RTIP=(R1(NR)+R1(NR-1))*0.5
      RHUB=(R1(1)+R1(2))*0.5
      IF(R1(1).EQ.R1(2))RHUB=0.
      DELR=RTIP-RHUB
      ISKIP = 0
      IF(IMOD.LT.4) ISKIP = 1
      DO 40 L=LS,LSS
      WRITE(6,300)J,L,(K,K=1,NTH)
      SP=(R1(L)-RHUB)/DELR
      DO 10 I=1,5
      WRITE(6,301)MSG(I,IMOD),(U(I,K,L),K=1,NTH)
      AV1=0.
      DO 5 K=1,NTH
      AV1=AV1+U(I,K,L)
5     CONTINUE
      AV(I)=AV1/NTH
10    CONTINUE
      WRITE(6,303)(MSG(I,IMOD),I=1,5),SP,(AV(I),I=1,5)
40    CONTINUE
300   FORMAT(///I4,I3,3X,17(5X,I2))
301   FORMAT(5X,A4,1X,17F7.3)
303   FORMAT(/,T2,'THETA AVERAGE PROPERTIES.',/,
     1T6,'%SPAN',5A7/
     2T5,6(F6.3,1X))
      IF(ISKIP.NE.0) GO TO 100
      WRITE(6,305)
305   FORMAT(/,T3,'L  RADIUS MACH NUMBERS         P  ON  T  ON',
     1'  ANGLES             ',/,T13,
     2'MERID  TANGEN TOTAL  PREF   TREF   WHIRL  MERID')
306   FORMAT(T2,I2,2(2X,F5.3),2X,F5.2,3(2X,F5.3),2(2X,F5.1),2X,F5.3)
C
      DO 90 L=1,NR
      R(L)=R1(L)
      DO 50 I=1,5
      AV1=0.
      DO 45 K=1,NTH
45    AV1=AV1+U(I,K,L)
      AV(I)=AV1/NTH
50    CONTINUE
      TMR=(AV(2)**2 + AV(4)**2)
      TMR=SQRT(TMR)
      TM=TMR**2 +AV(3)**2
      TM=SQRT(TM)
      AV(1)=AV(1)*GAMMA
      PTS(L)=AV(1)
      TTS(L)=AV(5)
      AV(5)=AV(5)-1.0
      ATW=ATAN(AV(3)/TMR)*TODEG
      ATR=ATAN(AV(2)/AV(4))*TODEG
      WRITE(6,306)L,R(L),TMR,AV(3),TM,AV(1),TTS(L),ATW,ATR
90    CONTINUE
C
      AV1=0.
      AV2=0.
      PTS(1)=PTS(1)*0.5
      TTS(1)=TTS(1)*0.5
      PTS(NR)=PTS(NR)*0.5
      TTS(NR)=TTS(NR)*0.5
      DO 95 L=2,NR1
      LP1=L+1
      LM1=L-1
      RU=(R(LP1)+R(L))*0.5
      RL=(R(L)+R(LM1))*0.5
      R2=RU**2-RL**2
      AV1=AV1+PTS(L)*R2
      AV2=AV2+TTS(L)*R2
95    CONTINUE
      AREA=(RTIP**2-RHUB**2)
      AV1=AV1/AREA
      AV2=AV2/AREA
C
100   WRITE(6,302)DMD
302   FORMAT(/,T2,'AREA INTEGRATED PROPERTIES.',/,
     1         T2,'TOTAL MASS  FLOW  = ',F6.2)
      IF(ISKIP.NE.0) GO TO 110
      WRITE(6,307)AV1,AV2
307   FORMAT(1X,T2,'TOTAL PRES. RATIO = ',F6.4,/,
     1          T2,'TOTAL TEMP. RATIO = ',F6.4)
110   CONTINUE
      RETURN
      END
      SUBROUTINE PROJR(U,K,L,INR,INZ,IPR,IPZ,VP,VN)
C
C              THIS ROUTINE FINDS THE NORMAL AND PARALLEL VELOCITIES
C
      REAL U(5,17,18),INR,INZ,IPR,IPZ
C
      VN=U(2,K,L)*INR +U(4,K,L)*INZ
      VN=VN/U(1,K,L)
C
C              FIND PARALLEL COMPONENTS
C
      IPZ=ABS(INR)
      IPR=-INZ*IPZ/INR
C
      VP=U(2,K,L)*IPR + U(4,K,L)*IPZ
      VP=VP/U(1,K,L)
C
      RETURN
      END
      SUBROUTINE PROJT(V,U,G1,G2,G3)
      COMMON/PASS/K,L,R
      COMMON/PROP/GAMMA,W,TP
      REAL U(5,17,18)
C
      RHOR= U(1,K,L)
      V= G1*U(2,K,L) + G3*U(4,K,L)
      V = V + G2*(U(3,K,L) - W*RHOR*R)
      V= V/RHOR
C
      RETURN
      END
      SUBROUTINE PTASS(IPRT,NM,JSS1,JSS2,LQ1,LQ2,ISEC,KSAV)
C
      COMMON/CONTRL/NSTEP,NOR(10,3)
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/PROP/GAMMA,W,G1
      COMMON/ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      COMMON/PASS/K,L,R
      COMMON/PROUT/PRPS(100),PRSS(100)
      REAL UA(120),U(5,17,18),R1(18),MDOT
      REAL TU1(1728)
      EQUIVALENCE (TU1,U)
      REAL XR1(18),XX1(18),RX1(18),RR1(18)
      REAL JAC(50)
      REAL MDSAV
      COMMON /MASS/MDSAV(100)
      EQUIVALENCE(UA(1),R1(1))
      EQUIVALENCE (UA(73),RR1(1))
      EQUIVALENCE (UA(19),XX1(1)),(UA(37),XR1(1)),
     *       (UA(55),RX1(1))
      REAL MSG(5,5),MX,M,SYS(5)
      DATA MSG/4HRHOP,4HURP ,4HUTHP,4HUXP ,4HETP ,
     *   4HRHO ,4HUR  ,4HUTH ,4HUX  ,4HET  ,
     *   4HRHO ,4HMR  ,4HMTH ,4HMX  ,4HP   ,
     *   4HPT  ,4HMR  ,4HMTH ,4HMX  ,4HTT  ,
     *   4HPT  ,4HMR  ,4HMTH ,4HMX  ,4HTT  /
      DATA SYS/4HCONS,4HABS ,4HREL ,4HREL ,4HABS /
C
C              MODE=1  =CONSERVATION VARIABLES
C              MODE=2  = PHYSICAL VARIABLES IN ABS SYS.
C              MODE=3  = RHO,MR,MTH,MX,P    IN REL SYSTEM
C              MODE=4  =PT,MR,MTK,MX,TT IN REL SYS.
C              MODE=5  =PT,MR,MTH,MX,TT IN ABS SYSTEM
C
      CALL IO(1,UA,0,1,TU1,JAC)
      JSTRT=JSS1
      JSTOP=JSS2
      LS=LQ1
      LSS=LQ2
      NY=NTH
      G2=G1/2.
      G3=GAMMA/G1
      MODE = IABS(NM)
C
      IF(NM.LT.0) GO TO 55
      WRITE(6,301)JSTRT,JSTOP,SYS(MODE)
301   FORMAT(1H1,T2,'GRID OUTPUT FOR',I4,' TO',I4,' IN ',A4,' SYSTEM')
      WRITE(6,304)(TITLE(I),I=1,80)
304   FORMAT(1X,T1,80A1)
55    IMOD = MODE
      IF(MODE.EQ.5)IMOD=5
      DO 500 J=JSTRT,JSTOP
      JJ=J
      IA=0
      CALL IO(1,UA,IA,JJ,TU1,JAC)
      JKK=3
      IF(J.EQ.1 .OR. J.EQ.NX)JKK=2
      IF(IMOD.EQ.1)JKK=0
      CALL MD(MDOT,U,R1,RR1,RX1,XX1,XR1,JAC,JKK)
C.... IF NM.LT.0, SAVE MASS FLOWS AND RETURN.
      MDSAV(J)=MDOT
      L=NR/2
      R=R1(L)
      K=2
      CALL FINDP(U,P,V)
      PRPS(J)=P
      K=NTH-1
      CALL FINDP(U,P,V)
      PRSS(J)=P
      IF(NM.LT.0) GO TO 500
      IF(IMOD.EQ.1)GO TO 400
      IF(IMOD.GT.2)GO TO 200
C
      DO 60 L=1,NR
      DO 60 K=1,NTH
      RHOR= U(1,K,L)
      U(1,K,L)= RHOR/R1(L)
      DO 60 I=2,5
60    U(I,K,L)=U(I,K,L)/RHOR
      GO TO 400
C
200   DO 90 L=1,NR
      R=R1(L)
      DO 90 K=1,NTH
      CALL FINDP(U,P,V)
      E=(U(5,K,L)-0.5*V)/U(1,K,L)
      T=GAMMA*G1*E
      IF(T.GT.0.)GO TO 936
      WRITE(6,510)J,K,L
510   FORMAT(' PRESSURE LESS THAN ZERO AT ',3I5)
      T=1.0
936   CONTINUE
      C=SQRT(T)
      RHOR=U(1,K,L)
      U(4,K,L)=U(4,K,L)/C/RHOR
      U(2,K,L)=U(2,K,L)/C/RHOR
      U(3,K,L)=U(3,K,L)/C/RHOR
      U(5,K,L)=P
      U(1,K,L)=RHOR/R
      IF(IMOD.EQ.3 .OR. IMOD.EQ.4)U(3,K,L)=U(3,K,L)-W*R/C
90    CONTINUE
      IF(IMOD.EQ.3) GO TO 400
      DO 100 L=1,NR
      DO 100 K=1,NTH
      M=0.
      DO 91 I=2,4
      M=M+U(I,K,L)**2
91    CONTINUE
      T1=1.0+G2*M
      PT=U(5,K,L)*T1**G3
      E=U(5,K,L)/U(1,K,L)/G1
      T=E*GAMMA*G1
      TT=T*T1
      U(5,K,L)=TT
      U(1,K,L)=PT
100   CONTINUE
400   CONTINUE
      CALL PRINT(U,IMOD,MSG,J,LS,LSS,MDOT,R1)
500   CONTINUE
      RETURN
      END
      SUBROUTINE SPLINT(N, MAX, X, Y, Z, YINT, DYDX, D2YDX, S, EM, WORK)
C
C        CALLING SEQUENCE
C              N - THE NUMBER OF SPLINE POINTS (MINIMUM 2)
C              MAX - THE NUMBER OF INTERPOLATION POINTS
C              X - VECTOR OF X COORDINATES OF N SPLINE POINTS
C              Y - VECTOR OF Y COORDINATES OF N SPLINE POINTS
C              Z - VECTOR OF X-COORDINATES OF INTERPOLATION POINTS
C              YINT - VECTOR OF Y COORDINATES CORRESPONDING TO POINTS Z
C              DYDX - VECTOR OF DY/DX         CORRESPONDING TO POINTS Z
C              D2YDX- VECTOR OF D2Y/DX2       CORRESPONDING TO POINTS Z
C              S,EM - ARRAYS OF DIMENSION N CALCULATES BY SPLINT
C              WORK - ARRAY USED BY SPLINT AS WORKING SPACE
C
      DIMENSION X(N), Y(N)
      DIMENSION Z(MAX), YINT(MAX), DYDX(MAX), D2YDX(MAX)
      DIMENSION S(N), EM(N)
      DIMENSION WORK(2,N)
1000  FORMAT(17H0OUT OF RANGE X =, G14.6)
1001  FORMAT('0   ERROR: SPLINT ENTERED WITH N LESS THAN TWO')
C***********************************************************************
      IF(N.LT.2) GO TO 900
      DO 10 I=2,N
10    S(I)=X(I)-X(I-1)
      NO=N-1
      WORK(1,1) =-.5
      WORK(2,1)=0.
      IF(2.GT.NO) GO TO 35
      DO 30 I=2,NO
      SI=S(I)
      SI1=S(I+1)
      F=(Y(I+1)-Y(I))/SI1 - (Y(I)-Y(I-1))/SI
      SID6=0.16666667*SI
      W=(SI+SI1)*.3333333-SID6*WORK(1,I-1)
      WORK(1,I)=(SI1*.1666667)/W
      WORK(2,I)=(F-SID6*WORK(2,I-1))/W
30    CONTINUE
35    EM(N)=(.5*WORK(2,NO))/(1.+.5*WORK(1,NO))
      DO 40 I=2,N
      K=N+1-I
40    EM(K)=WORK(2,K)-WORK(1,K)*EM(K+1)
      IF(MAX.LT.1) RETURN
      DO 150 I=1,MAX
      ZI=Z(I)
      K=2
      IF(ZI-X(1)) 50, 100, 55
50    IF(ZI.LT.(1.1*X(1)-.1*X(2))) WRITE(6,1000) ZI
      GO TO 100
55    K=N
      IF(ZI-X(N)) 65, 100, 60
60    IF(ZI.GT.(1.1*X(N)-.1*X(N-1))) WRITE(6,1000) ZI
      GO TO 100
65    MN1=2
      MX=N
70    K=(MX+MN1)/2
      IF(MX.EQ.MN1) GO TO 100
      IF(ZI-X(K)) 75, 100, 80
75    MX=K
      GO TO 70
80    MN1=K+1
      GO TO 70
100   XZ=X(K)-ZI
      XZ2=XZ*XZ
      ZX=ZI-X(K-1)
      ZX2=ZX*ZX
      SK=S(K)
      SK2=SK*2.
      SKD6=SK*.16666667
      EMK=EM(K)
      EMK1=EM(K-1)
      YKDS=Y(K)/SK
      YK1DS=Y(K-1)/SK
      YINT(I)=(EMK1*XZ2*XZ+EMK*ZX2*ZX)*.1666667/SK + (YKDS-EMK*SKD6)*ZX
     *   + (YK1DS-EMK1*SKD6)*XZ
      DYDX(I)= (EMK*ZX2-EMK1*XZ2)/SK2 + YKDS - YK1DS - (EMK-EMK1)*SKD6
      D2YDX(I)=(EMK1*XZ+EMK*ZX)/SK
150   CONTINUE
      RETURN
900   WRITE(6,1001)
      RETURN
      END
      SUBROUTINE  START(IBEG)
      COMMON/UPBOND/PDOWN,UPMACH
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      REAL USAV(18)
      COMMON/PROP/GAMMA,W,G1
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/US1/U1(1728,100),UA1(120,100),JAC1(20,100)
      REAL U(5,17,18),TU(1728),UA(120),JAC(50),JAC1
      COMMON/TGSAV/TTSAV(18,18,100),TRSAV(18,18,100),
     * TXSAV(18,18,100),CURSAV(18,2,2,100),
     * COSNSV(18,3,2,100)
      EQUIVALENCE (U,TU)
      COMMON /FLIPS/IFLIPX,IFLIPT,IFLIPR,ICYCLE
      COMMON/DIST/TDA(20),RREF(20),UTHETA(20),INOPT
      HO=1./G1
      IF(IBEG.EQ.1) GO TO 15
C.... READ OLD SOLUTION.
      READ(1) DUM,DUM,DUM,DUM,DUM,DUM
      DO 10 J=1,NX
      READ(1)TU
      DO 5 I=1,1728
5     U1(I,J)=TU(I)
10    CONTINUE
      READ(1,END=11)IFLIPX,IFLIPT,IFLIPR,ICYCLE
11    CONTINUE
      REWIND 1
      GO TO 1000
C
C.... GENERATE STARTING SOLUTION.
15    CONTINUE
      RG1=1./G1
      RGAMG1=1./(GAMMA*G1)
      T=1.0+(UPMACH**2)*(G1)/2.
      T=1.0/T
      UX=UPMACH*SQRT(T)
      WRITE(6,600) UX
600   FORMAT(T2,'UX = ',F10.5)
      DO 100 J=1,NX
      DO 60 L=1,NR
      IF(J.GE.N1)GO TO 20
C     UPSTREAM
      USAV(L)=0.0
      GO TO 30
20    IF(J.GT.N2-2)GO TO 30
C     INSIDE BLADE ROW OR DOWNSTREAM
      BETA=TXSAV(2,L,J)/TTSAV(2,L,J)
      USAV(L)=(W-BETA*UX)*UA1(L,J)
30    UTH=UTHETA(L)+USAV(L)
      R=UA1(L,J)
      U2=(UX*UX+UTH*UTH)*0.5
      H=HO+W*R*UTH-U2
      T=G1*H
      RHO=T**RG1
      E=T*RGAMG1
      ET=E+U2
      RHOR=RHO*R
      DO 55 K=1,NTH
      U(1,K,L)=RHOR
      U(2,K,L)=0.0
      U(3,K,L)=UTH*RHOR
      U(4,K,L)=UX*RHOR
      U(5,K,L)=ET*RHOR
55    CONTINUE
60    CONTINUE
      DO 65 I=1,1728
65    U1(I,J)=TU(I)
100   CONTINUE
      IFLIPX=1
      IFLIPT=1
      IFLIPR=1
      ICYCLE=0
      CALL TPM
C
1000  CONTINUE
      RETURN
      END
      SUBROUTINE STEPR(U2,R1,RR1,RX1,IFW,INX,UT,RK)
C
C     THIS ROUTINE DOES THE R-TIME STEP INTEGRATION
C     THIS VERSION DOES FORWARD-BACKWARD FLIPPING
C
      REAL UT(5,17,18),U2(5,17,18)
      REAL F1(5),F2(5),F3(5),G1(5),G2(5),G3(5)
      REAL H1(5),H2(5),H3(5),K1(5),K2(5),K3(5)
      REAL R1(18),RR1(18),RX1(18)
      REAL VISP(18),VIS(5,18)
      REAL JAC1,RK(5)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /VISCOS/CAP,KA(5)
      COMMON /PASS/K,L,R
      COMMON /DPM/NRLOW,NRHIGH
      COMMON /FLIP/IFLIP
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXTPE
      COMMON /RESIDS/RESID(18)
      COMMON /US1/UTEM(1728,100),UA1(120,100),JAC1(20,100)
      COMMON /ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
      NR1  =NRHIGH-1
      NRL  =NRLOW+1
      NR1FP=NR1*IFLIP
      NR1FM=-NR1FP
      LLPFP=NRL*IFLIP
      LLPFM=-LLPFP
      NTH1 =NTH-1
      TAURF=TAUR*IFLIP
C
      DO 3 K=1,NTH
      DO 3 I=1,5
      UT(I,K,NRLOW )=U2(I,K,NRLOW)
3     UT(I,K,NRHIGH)=U2(I,K,NRHIGH)
C
      DO 30 K=2,NTH1
C
C     VISP=STATIC PRESSURE FOR ARTIFICIAL VISCOSITY
C
      L=NRHIGH
      R=R1(L)
      CALL FINDP(U2,VISP(L),VV)
C
      L=NRLOW
      R=R1(L)
      CALL FINDP(U2,VISP(L),VV)
C
C     FOR IFLIP=+1
C     PREDICTOR : UT = U - DT*DH/DR-DT*DF/DR
C     OPPOSITE FOR IFLIP=-1
C
      DO 10 LL=NRL,NR1
      LP   =LL+IFLIP
      RR1LL=RR1(LL)
      RR1LP=RR1(LP)
      RX1LL=RX1(LL)
      RX1LP=RX1(LP)
      DIP  =JAC1(LL,INX)/JAC1(LP,INX)
C
      L=LL
      R=R1(LL)
      CALL CAFGHK(U2,F2,G2,H2,K2,VISP(LL),V)
C
      L=LP
      R=R1(LP)
      CALL CAFGHK(U2,F3,G3,H3,K3,PP3,V)
C
      IF(LL.NE.NR1FP.AND.LL.NE.LLPFM)GO TO 6
      H0=H2(5)/H2(1)
      H3(5)=H0*H3(1)
      F3(5)=H0*F3(1)
6     CONTINUE
C
      IF(IFCL.NE.1) GO TO 1007
      DO 7 I=1,5
      CAPF2= F2(I)*RR1LL+H2(I)*RX1LL
      CAPF3=(F3(I)*RR1LP+H3(I)*RX1LP)*DIP
      T=CAPF3-CAPF2
7     UT(I,K,LL)=U2(I,K,LL)-TAURF*T+DT*RK(I)*K2(I)
      GO TO 10
1007  CONTINUE
      DO 8 I=1,5
      T=RX1LL*(H3(I)-H2(I))+RR1LL*(F3(I)-F2(I))
8     UT(I,K,LL)=U2(I,K,LL)-TAURF*T+DT*RK(I)*K2(I)
10    CONTINUE
C
C     VIS = ARTIFICIAL VISCOSITY
C
      DO 15 LL=NRL,NR1
      LLP=LL+1
      LLM=LL-1
      RJP=R1(LL)/R1(LLP)
      RJM=R1(LL)/R1(LLM)
      DJ=ABS(VISP(LLP)-2.*VISP(LL)+VISP(LLM))
     1     /(VISP(LLP)+2.*VISP(LL)+VISP(LLM))
      DO 15 I=1,5
15    VIS(I,LL)=CAP*DJ*(RJP*U2(I,K,LLP)-2.*U2(I,K,LL)+RJM*U2(I,K,LLM))
C
C     FOR IFLIP=+1
C     CORRECTOR STEP : U(N+1)=0.5(U+UT-DT*(DH/DR+DF/DR)+VIS)
C     OPPOSITE FOR IFLIP=-1.  CALCULATION ALWAYS CENTERED AT 2.
C
      DO 20 LL=NRL,NR1
      LM   =LL-IFLIP
      RR1LM=RR1(LM)
      RR1LL=RR1(LL)
      RX1LM=RX1(LM)
      RX1LL=RX1(LL)
      DIM  =JAC1(LL,INX)/JAC1(LM,INX)
C
      L=LM
      R=R1(LM)
      CALL CAFGHK(UT,F1,G1,H1,K1,PP1,V)
C
      L=LL
      R=R1(LL)
      CALL CAFGHK(UT,F2,G2,H2,K2,PP,V)
C
      IF(LL.NE.NR1FM.AND.LL.NE.LLPFP)GO TO 16
      H0=H2(5)/H2(1)
      F1(5)=H0*F1(1)
      H1(5)=H0*H1(1)
      UT(5,K,1)=UT(5,K,3)/R1(3)*R1(1)
16    CONTINUE
C
      IF(IFCL.NE.1) GO TO 1017
      DO 17 I=1,5
      CAPF1=(F1(I)*RR1LM+H1(I)*RX1LM)*DIM
      CAPF2= F2(I)*RR1LL+H2(I)*RX1LL
      T=CAPF2-CAPF1
      UCOR=U2(I,K,LL)-TAURF*T+DT*RK(I)*K2(I)
17    U2(I,K,LL)=.5*(UT(I,K,LL)+UCOR+VIS(I,LL))
      GO TO 20
1017  CONTINUE
      DO 18 I=1,5
      T=RX1LL*(H2(I)-H1(I))+RR1LL*(F2(I)-F1(I))
      UCOR=U2(I,K,LL)-TAURF*T+DT*RK(I)*K2(I)
18    U2(I,K,LL)=.5*(UT(I,K,LL)+UCOR+VIS(I,LL))
20    CONTINUE
30    CONTINUE
C     RESID = T.E. STATIC PRESSURES
      IF(INX.NE.N2)RETURN
      K = NTH1
      DO 700 L=NRLOW,NRHIGH
      R = R1(L)
700   CALL FINDP(U2,RESID(L),V)
      RETURN
      END
      SUBROUTINE STEPT(U2,UT,INX,IFW,TT,TR,TX,R1,RK)
C
C     THIS SUBRUOTINE DOES BOTH THE FORWARD AND BACKWARD THETA
C     STEP IN ONE CALL.
C     FORWARD-BAKCWARD FLIPPING IS DONE
C
      REAL FSAV(5,4,17)
      REAL U2(5,17,18),UT(5,17,18),TT(16),TR(15,16),TX(15,16),R1(18)
      REAL F2(5),G2(5),H2(5),K2(5)
      REAL VISP(18),VIS(5,18),RK(5)
      REAL TRF(17),TXF(17),GHAT(5,17)
C
      COMMON/PASS/K,L,R
      COMMON/VISCOS/CAP,KA(5)
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/FLIP/IFLIP
      COMMON/ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
C
      NR1 =NR-1
      NTH1=NTH-1
      TAUTF=TAUT*IFLIP
C
C     DO FORWARD AND BACKWARD STEP IN ORDER SPECIFIED BY IFLIP.
C     IFLIP = 1 => FORWARD, BACKWARD
C     IFLIP =-1 => BACKWARD, FORWARD
C     CALCULATION IS ALWAYS CENTERED ON 2.  DEPENDING ON IFLIP
C     3 IS AHEAD OR BEHIND 2.
C
      DO 120 LQQ=2,NR1
      L  =LQQ
      L1 =L-1
      TTL=TT(L1)
      TTR=TTL
      R  =R1(L)
C
C     THETA METRICS FOR CONSERVATIVE FORM OF EQUATIONS
C
      IF(IFCL.NE.1) GO TO 2
      DO 1 K=2,NTH1
      KM=K-1
      TRF(K)=TR(KM,L1)
1     TXF(K)=TX(KM,L1)
      TXT     =TXF(3)-TXF(2)
      TXF(1)  =TXF(2)   -TXT
      TXF(NTH)=TXF(NTH1)+TXT
      TRT     =TRF(3)-TRF(2)
      TRF(1)  =TRF(2)   -TRT
      TRF(NTH)=TRF(NTH1)+TRT
2     CONTINUE
C
      DO 4 K=1,NTH
      CALL CAFGHK(U2,F2,G2,H2,K2,PP,V)
      VISP(K)=PP
      DO 4 I=1,5
      FSAV(I,1,K)=F2(I)
      FSAV(I,2,K)=G2(I)
      FSAV(I,3,K)=H2(I)
4     FSAV(I,4,K)=K2(I)
      IF(IFCL.NE.1)GO TO 1010
      DO 1009 K=1,NTH
      DO 1009 I=1,5
      GHAT(I,K)=FSAV(I,1,K)*TRF(K)+FSAV(I,2,K)*TTR
     1         +FSAV(I,3,K)*TXF(K)
1009  CONTINUE
1010  CONTINUE
C
C     ARTIFICIAL VISCOSITIES
C
      DO 5 K=2,NTH1
      KP=K+1
      KM=K-1
      V1=ABS(VISP(KP)-2.*VISP(K)+VISP(KM))
      V2=    VISP(KP)+2.*VISP(K)+VISP(KM)
      DJ=CAP*V1/V2
      DO 5 I=1,5
5     VIS(I,K)=DJ*(U2(I,KP,L)-2.*U2(I,K,L)+U2(I,KM,L))
C
C     PREDICTOR STEP
C
      DO 10 KK=2,NTH1
      KPLUS=KK+IFLIP
      K    =KK
C
      IF(IFCL.NE.1)GO TO 1006
      DO 6 I=1,5
      T=GHAT(I,KPLUS)-GHAT(I,K)
6     UT(I,K,L) = U2(I,K,L)-TAUTF*T+DT*RK(I)*FSAV(I,4,K)
      GO TO 10
C
1006  CONTINUE
      KM   =K-1
      TXL  =TX(KM,L1)
      TRL  =TR(KM,L1)
      DO 7 I=1,5
      T=TTL*(FSAV(I,2,KPLUS)-FSAV(I,2,K))
     1 +TXL*(FSAV(I,3,KPLUS)-FSAV(I,3,K))
     2 +TRL*(FSAV(I,1,KPLUS)-FSAV(I,1,K))
7     UT(I,K,L) = U2(I,K,L)-TAUTF*T+DT*RK(I)*FSAV(I,4,K)
10    CONTINUE
C
      DO 1004 K=1,NTH
      CALL CAFGHK(UT,F2,G2,H2,K2,PP,V)
      DO 1004 I=1,5
      FSAV(I,1,K)=F2(I)
      FSAV(I,2,K)=G2(I)
      FSAV(I,3,K)=H2(I)
1004  FSAV(I,4,K)=K2(I)
      IF(IFCL.NE.1)GO TO 1020
      DO 1019 K=1,NTH
      DO 1019 I=1,5
      GHAT(I,K)=FSAV(I,1,K)*TRF(K)+FSAV(I,2,K)*TTR
     1         +FSAV(I,3,K)*TXF(K)
1019  CONTINUE
1020  CONTINUE
C
C     CORRECTOR STEP
C
      DO 110 KK=2,NTH1
      KMIN=KK-IFLIP
      K   =KK
C
      IF(IFCL.NE.1)GO TO 1012
      DO 12 I=1,5
      T=GHAT(I,K)-GHAT(I,KMIN)
      UCOR=U2(I,K,L)-TAUTF*T+DT*RK(I)*FSAV(I,4,K)
12    U2(I,K,L)=.5*(UT(I,K,L)+UCOR+VIS(I,K))
      GO TO 110
C
1012  CONTINUE
      KM  =KK-1
      TXL=TX(KM,L1)
      TRL=TR(KM,L1)
      DO 13 I=1,5
      T=TTL*(FSAV(I,2,K)-FSAV(I,2,KMIN))
     1 +TXL*(FSAV(I,3,K)-FSAV(I,3,KMIN))
     2 +TRL*(FSAV(I,1,K)-FSAV(I,1,KMIN))
      UCOR=U2(I,K,L)-TAUTF*T+DT*RK(I)*FSAV(I,4,K)
13    U2(I,K,L)=.5*(UT(I,K,L)+UCOR+VIS(I,K))
110   CONTINUE
120   CONTINUE
C
      RETURN
      END
      SUBROUTINE STEPX(U1,U2,U3,UP1,R1,R2,R3,XR,XX,IFW,INX)
C
C     THIS ROUTINE DOES THE X-TIME STEP INTEGRATION
C     IFLIP=+1 => J=1,2,3...
C     IFLIP=-1 => J=100,99,98...
C     CALCULATION CENTERED ON J2
C     INPUT:  U1,U2,U3 = SOLUTIONS AT J1,J2,J3
C             UP1      = PREDICTED U AT J1 FROM LAST STEPX CALL
C     CODE:   UP2      = PREDICTED U AT J2 = F(U2,U3)
C             VIS      = ART. VISC.  AT J2 = F(U1,U2,U3)
C     OUTPUT: U1       = CORRECTED U AT J2 = F(UP1,UP2)
C             UP1      = UP2 = PREDICTED U AT J2
C
      REAL UP1(5,17,18),UP2(5)
      REAL U1(5,17,18),U2(5,17,18),U3(5,17,18)
      REAL R1(18),R2(18),R3(18)
      REAL XR(18),XX(18)
      REAL F1(5),F2(5),F3(5),H1(5),H2(5),H3(5)
      REAL K1(5),K2(5),K3(5),G1(5),G2(5),G3(5)
      REAL VIS(5),JAC1
C
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /VISCOS/CAP,KA(5)
      COMMON /PASS/K,L,R
      COMMON /FLIP/IFLIP
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXTPE
      COMMON /ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
      COMMON /US1/UTEM(1728,100),UA1(120,100),JAC1(20,100)
C
      NR1=NR-1
      NT1=NTH-1
      INXP=INX+IFLIP
      INXM=INX-IFLIP
      TAUXF=TAUX*IFLIP
C
C     DO BOTH FORWARD AND BACKWARD STEP AS SPECIFIED BY IFLIP
C
      DO 100 L=2,NR1
      XRL=XR(L)
      XXL=XX(L)
      LXX=18+L
      LXR=36+L
      XR3=UA1(LXR,INXP)
      XR1=UA1(LXR,INXM)
      XX3=UA1(LXX,INXP)
      XX1=UA1(LXX,INXM)
      D1=JAC1(L,INX)/JAC1(L,INXM)
      D3=JAC1(L,INX)/JAC1(L,INXP)
C
      DO 100 K=2,NT1
C
C     PREDICTOR: UT = U - DT*DH/DZ  -DT*DF/DZ
C
      R=R1(L)
      CALL FINDP(U1,PP1,VV)
C
      R=R2(L)
      CALL CAFGHK(U2,F2,G2,H2,K2,PP2,V)
C
      R=R3(L)
      CALL CAFGHK(U3,F3,G3,H3,K3,PP3,V)
      RJP=R2(L)/R3(L)
C
      IF(IFCL.NE.1) GO TO 1007
      DO 7 I=1,5
      CAPH2= F2(I)*XRL+H2(I)*XXL
      CAPH3=(F3(I)*XR3+H3(I)*XX3)*D3
7     UP2(I)=U2(I,K,L)-TAUXF*(CAPH3-CAPH2)
      GO TO 10
1007  CONTINUE
      DO 8 I=1,5
      T=XXL*(H3(I)-H2(I))+XRL*(F3(I)-F2(I))
8     UP2(I)=U2(I,K,L)-TAUXF*T
10    CONTINUE
C
C     CORRECTOR: U(N+1)=0.5(U+UT-DT*(DH/DZ+DF/DZ)+VISC)
C
      R=R1(L)
      CALL CAFGHK(UP1,F1,G1,H1,K1,PP,V)
C
      DO 20 I=1,5
20    UP1(I,K,L)=UP2(I)
C
      R=R2(L)
      CALL CAFGHK(UP1,F2,G2,H2,K2,PP,V)
      RJM=R2(L)/R1(L)
C
C     ARTIFICIAL VISCOSITY TERMS
      A1=ABS(PP3-2.*PP2+PP1)
      A2=    PP3+2.*PP2+PP1
      AJ=CAP*A1/A2
      DO 12 I=1,5
      AJP=RJP*U3(I,K,L)-U2(I,K,L)
      AJM=U2(I,K,L)-RJM*U1(I,K,L)
12    VIS(I)=AJ*(AJP-AJM)
C
      IF(IFCL.NE.1) GO TO 1013
      DO 13 I=1,5
      CAPH1=(F1(I)*XR1+H1(I)*XX1)*D1
      CAPH2= F2(I)*XRL+H2(I)*XXL
      T=CAPH2-CAPH1
13    U1(I,K,L)=0.5*(U2(I,K,L)+UP2(I)-TAUXF*T+VIS(I))
      GO TO 100
1013  CONTINUE
      DO 14 I=1,5
      T=XXL*(H2(I)-H1(I))+XRL*(F2(I)-F1(I))
14    U1(I,K,L)=0.5*(U2(I,K,L)+UP2(I)-TAUXF*T+VIS(I))
C
100   CONTINUE
C
C     COPY BOUNDARY VALUES INTO NEW MATRIX
C
      DO 200 L=2,NR1
      DO 200 I=1,5
      U1(I,1,L)=U2(I,1,L)
200   U1(I,NTH,L)=U2(I,NTH,L)
C
      DO 300 K=2,NT1
      DO 300 I=1,5
      U1(I,K,1)=U2(I,K,1)
300   U1(I,K,NR)=U2(I,K,NR)
C
      INX=INX+IFLIP
      RETURN
      END
      SUBROUTINE TPM
      COMMON/CONTRL/NSTEP,NOR(10,3)
      COMMON/TQ/UA(120),TU(1728),JAC(50)
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      REAL JAC
C
      IR=1
      XM=100
      TM=100
      RM=100
      DO 100 J=1,NX
      CALL IO(IR,UA,0,J,TU,JAC)
      CALL TSTIM(X,T,R,J)
      IF(X.GT.XM) GO TO 10
      XM=X
      JX=J
10    IF(T.GT.TM) GO TO 20
      TM=T
      JT=J
20    IF(R.GT.RM) GO TO 30
      RM=R
      JR=J
30    CONTINUE
100   CONTINUE
C
      WRITE(6,610)
610   FORMAT(//,'  CFL STABILITY LIMIT ON TIME STEP.')
      WRITE(6,600)XM,JX,TM,JT,RM,JR
600   FORMAT('  MIN X STEP OF ',F10.5,' AT J = ',I5,/,
     1       '  MIN T STEP OF ',F10.5,' AT J = ',I5,/,
     2       '  MIN R STEP OF ',F10.5,' AT J = ',I5)
      DTMAX = AMIN1(XM,TM,RM)
      CFM   = 0.95
      DTSUG = CFM*DTMAX
      WRITE(6,620)DTMAX,DTSUG
620   FORMAT('  DT MUST BE LESS THAN ',F7.5,/
     1       '  RECOMMENDED DT IS    ',F7.5)
C
      IF(DT.GT.0.0) GO TO 200
      DT=DTSUG
      DO 150 I=1,10
      DO 140 J=2,3
140   NOR(I,J)=1
150   NOR(I,1)=0
      NOR(3,2)=2
      NOR(1,1)=2
      NOR(2,1)=1
      NOR(3,1)=3
      NOR(4,1)=1
      NOR(5,1)=2
      TAUX=DT/DX
      TAUR=DT/DR
      TAUT=DT/DTH
C
200   RETURN
      END
      SUBROUTINE TRFIL(TT,TR,TX,CUR,COSN,J)
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/TGSAV/TTSAV(18,18,100),TRSAV(18,18,100),TXSAV(18,18,100),
     *CURSAV(18,2,2,100),COSNSV(18,3,2,100)
      REAL TT(16),TR(15,16),TX(15,16)
      REAL CUR(16,2,2),COSN(16,3,2)
      NTH1=NTH-1
      NR1=NR-1
      DO  100 L=2,NR1
      L1=L-1
      TT(L1)=TTSAV(1,L,J)
      DO  10 K=2,NTH1
      K1=K-1
      TR(K1,L1)=TRSAV(K,L,J)
      TX(K1,L1)=TXSAV(K,L,J)
10    CONTINUE
      DO  20 KX=1,2
      DO  30 KY=1,3
30    COSN(L1,KY,KX)=COSNSV(L,KY,KX,J)
      DO  20 KY=1,2
20    CUR(L1,KY,KX)=CURSAV(L,KY,KX,J)
100   CONTINUE
      RETURN
      END
      SUBROUTINE TSTIM(XMIN,TMIN,RMIN,JX)
      COMMON /GRID/ NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      REAL UA(120),U(5,17,18),TU1(1728)
      REAL R1(18),XX(18),XR(18),RR(18),RX(18)
      EQUIVALENCE (TU1,U),(UA,R1),(UA(19),XX),(UA(37),XR),
     *(UA(55),RX),(UA(73),RR)
      COMMON/TGSAV/TTSAV(18,18,100),TRSAV(18,18,100),
     *TXSAV(18,18,100),CURSAV(18,2,2,100),COSNSV(18,3,2,100)
      COMMON /PROP/GAMMA,W,G1
      COMMON /PASS/K,L,R
      COMMON /TQ/UA,TU1
      NR1=NR-1
      XMIN=10.
      TMIN=10.
      RMIN=10.
      NTH1=NTH-1
      G2=GAMMA*G1
C
      J=JX
C
C
      DO 50 L=1,NR
      R=R1(L)
C
C
C
      DO 45 K=1,NTH
      CALL FINDP(U,P,V)
      RHOR=U(1,K,L)
      E=(U(5,K,L)-0.5*V)/RHOR
      C=SQRT(GAMMA*G1*E)
      UX=U(4,K,L)/RHOR/C
      UR=U(2,K,L)/RHOR/C
      UT=U(3,K,L)/RHOR
      UT=(UT-W*R)/C
      TR=TRSAV(K,L,J)
      TX=TXSAV(K,L,J)
      TT=TTSAV(K,L,J)
      TT=TT/R
      CX=SQRT(XX(L)**2 +XR(L)**2)
      CR=SQRT(RR(L)**2 +RX(L)**2)
      CT=SQRT(TX**2+TT**2+TR**2)
      TXX=XX(L)*UX + XR(L)*UR
      TXX=ABS(TXX)+CX
      TXX=DX/TXX
      IF(TXX.LT.XMIN) XMIN=TXX
      TRR=RR(L)*UR + RX(L)*UX
      TRR=ABS(TRR)+CR
      TRR=DR/TRR
      IF(TRR.LT.RMIN) RMIN=TRR
      TP=TT*UT + TX*UX + TR*UR
      TP=ABS(TP)+CT
      TP=DTH/TP
      IF(TP.LT.TMIN) TMIN=TP
45    CONTINUE
50    CONTINUE
      RETURN
      END
      SUBROUTINE UNPROJ(UR,UZ,VP,VN,INR,INZ,IPR,IPZ)
C
C              THIS ROUTINE FINDS THE AXIAL AND RADIAL
C              VELOCITIES FROM THE NORMAL AND PARALLEL ONES
C
      REAL INR,INZ,IPR,IPZ
C
C
      T=VP*INR-VN*IPR
      D=IPZ*INR-IPR*INZ
      UZ=T/D
      UR=(VN-INZ*UZ)/INR
      RETURN
      END
      SUBROUTINE UNPROT(V1,V2,V3,UR,UTH,UZ)
      COMMON/PASSCO/GR,GT,GZ,AR,AT,AZ,BR,BT,BZ
C
C       THIS ROUTINE FINDS THE VELOCITY COMPONIENTS
C
      DET=BT*GZ-GT*BZ-(AT*GZ*BR-GR*AT*BZ)/AR
      T1=GZ*V3-BZ*V1
      T2=(GZ*BR-GR*BZ)*V2/AR
      UTH=(T1-T2)/DET
      UR=(V2-AT*UTH)/AR
      UZ=(V3-BR*UR-BT*UTH)/BZ
      RETURN
      END
      SUBROUTINE WALLP(RA,U,I)
      REAL RA(18),U(5,17,18)
      COMMON/PASS/K,L,R
      COMMON /PRESS/P(3,18,4),LR,IL,IH
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      INTEGER LRA(4)
      LRA(1)=2
      LRA(2)=NR-1
      LRA(3)=IH
      LRA(4)=IL
      NTH1=NTH-1
      DO 100 K=2,NTH1
      DO 50 LR=1,4
      L=LRA(LR)
      R=RA(L)
      CALL FINDP(U,P1,V)
      P(I,K,LR)=P1
50    CONTINUE
100   CONTINUE
      RETURN
      END
C   GENERAL PLOTTING PACKAGE FOR PRINTER OR TERMINAL
C
      SUBROUTINE PPINIT(NTYPE,IFLU,XMN,XMX,YMN,YMX)
C
C   INITIALIZATION ROUTINE FOR PRINTER PLOTTING PACKAGE
C
C
C      WHERE:
C
C     NTYPE - PLOT TYPE: 1 = LINEAR
C                        2 = SCALED SEMI-LOG (Y LINEAR)
C                        3 = SCALED SEMI-LOG (X LINEAR)
C                        4 = SCALED LOG-LOG
C                        5 = SEMI-LOG (Y LINEAR)
C                        6 = SEMI-LOG (X LINEAR)
C                        7 = LOG-LOG
C      IFLU - FORTRAN LOGICAL UNIT USED FOR PLOT OUTPUT
C       XMN - MIN X TO ACCEPT IN PLOT
C       XMX - MAX X FOR PLOT
C       YMN - MIN Y TO ACCEPT IN PLOT
C       YMX - MAX Y FOR PLOT
C
C             IF XMN AND XMX OR YMN AND YMX ARE ZERO THEN THE
C             PLOT IS NOT FORCED, AND THE PACKAGE WILL FIND THE
C             RANGES FOR X AND/OR Y FROM THE STORED DATA.
C
C
      COMMON/PPBLK/ ITYPE,IWRITE,XMIN,XMAX,YMIN,YMAX,ICALL,IPOS,
     *              KJ,IFIN,IC(50),IND(50),INUM(50),FORCEX,FORCEY
      LOGICAL FORCEX,FORCEY
C
C
C   INITIALIZE
C
      IWRITE=IFLU
      ITYPE=NTYPE
      XMIN=XMN
      XMAX=XMX
      YMIN=YMN
      YMAX=YMX
      FORCEX=XMIN.NE.0. .OR. XMAX.NE.0.
      FORCEY=YMIN.NE.0. .OR. YMAX.NE.0.
      IF(.NOT. FORCEY) GO TO 1
      IF(ITYPE.EQ.1.OR.ITYPE.EQ.2.OR.ITYPE.EQ.5) GO TO 1
      IF(YMIN.GT.0. .AND. YMAX.GT.0.) GO TO 1
      WRITE(IWRITE,100)
      FORCEY=.FALSE.
1     IF(.NOT. FORCEX) GO TO 3
      IF(ITYPE.EQ.1.OR.ITYPE.EQ.3.OR.ITYPE.EQ.6) GO TO 3
      IF(XMIN.GT.0. .AND. XMAX.GT.0.) GO TO 3
      FORCEX=.FALSE.
      WRITE(IWRITE,100)
100   FORMAT(/' **WARNING** FORCE PARAMETERS .LE. 0, PLOT NOT FORCED')
3     ICALL=0
      IPOS=1
      KJ=1
      IFIN=-1
      RETURN
      END
      SUBROUTINE PPLOT(XLEN,PHITE,ITITX,ITITY,WORK)
C
C   ROUTINE THAT PERFORMS THE PLOTTING FOR THE PACKAGE
C
C
C      WHERE:
C
C      XLEN - THE PLOT LENGTH IN INCHES (MAX - 11)
C     PHITE - THE PLOT HEIGHT IN INCHES
C     ITITX - A VECTOR CONTAINING X AXIS ANNOTATION IN A1 FORMAT.
C             A ZERO (0) INDICATES END OF TITLE.
C     ITITY - A VECTOR CONTAINING Y AXIS ANNOTATION IN A1 FORMAT.
C             A ZERO (0) INDICATES END OF TITLE.
C      WORK - THE SAME WORK VECTOR USED IN THE PPSTOR STAGE.
C
C
      DIMENSION ITITX(1),ITITY(1),WORK(1),XLABEL(12)
      DIMENSION LINE(112),CODES(6)
      DIMENSION PLC(11),FORMAT(7)
      COMMON/PPBLK/ ITYPE,IWRITE,XMIN,XMAX,YMIN,YMAX,ICALL,IPOS,
     *              KJ,IFIN,IC(50),IND(50),INUM(50),FORCEX,FORCEY
      LOGICAL FORCEX,FORCEY
C
      DATA FORMAT/'(1X,','A1,F','XXXX',',1X,','1H+,','XXXX',
     + 'A1) '/
      DATA CODES/'6.5','6.4','6.3','6.2','6.1','6.0'/
      DATA PLC/'  12','  22','  32','  42','  52'
     * ,'  62','  72','  82','  92',' 102',' 112'/
      DATA MIX/'X'/,IMNS/'-'/,IOH/'O'/,IBLK/' '/,IPLS/'+'/
C  CLIP PLOT AT SINGLE PRECISION MAXIMUM
      DATA SPMAX/1.E38/
C
C  PLOT STORED VALUES
C
      IF(IPOS .GT. 1) GO TO 8
      WRITE(IWRITE,702)
702   FORMAT(/' ***ERROR*** NO PLOTS STORED')
      RETURN
C
C   GET X AND Y AXES SET UP
C
8     IF(.NOT. FORCEX) GO TO 81
      IF(ITYPE.NE.2.AND.ITYPE.NE.4.AND.ITYPE.NE.5.AND.ITYPE.NE.7)
     * GO TO 81
      XMIN=ALOG10(XMIN)
      XMAX=ALOG10(XMAX)
81    IF(.NOT. FORCEY) GO TO 80
      IF(ITYPE.NE.3.AND.ITYPE.NE.4.AND.ITYPE.NE.6.AND.ITYPE.NE.7)
     * GO TO 80
      YMIN=ALOG10(YMIN)
      YMAX=ALOG10(YMAX)
C   FIND X AND Y RANGE
80    IF(FORCEX .AND. FORCEY) GO TO 82
      IF(.NOT. FORCEY) YMAX=-SPMAX
      IF(.NOT. FORCEY) YMIN=SPMAX
      IF(.NOT. FORCEX) XMIN=SPMAX
      IF(.NOT. FORCEX) XMAX=-SPMAX
      DO 83 I=1,IFIN,2
      IF(.NOT. FORCEY) YMAX=AMAX1(YMAX,WORK(I+1))
      IF(.NOT. FORCEY) YMIN=AMIN1(YMIN,WORK(I+1))
      IF(.NOT. FORCEX) XMAX=AMAX1(XMAX,WORK(I))
      IF(.NOT. FORCEX) XMIN=AMIN1(XMIN,WORK(I))
83    CONTINUE
C
82    ILEN=XLEN+.49
      ILEN=MAX0(ILEN,1)
      ILEN=MIN0(ILEN,11)
      PLEN=ILEN*10
      XMX=XMAX
      XMN=XMIN
      YMX=YMAX
      YMN=YMIN
      IF(ITYPE.NE.3.AND.ITYPE.NE.4) GO TO 74
      IF(YMN .GE. 0.) YMIN=IFIX(YMN)
      IF(YMN .LT. 0.) YMIN=IFIX(YMN)-1
      IF(YMX .GE. 0.) YMAX=IFIX(YMX)+1
      IF(YMX .LT. 0.) YMAX=IFIX(YMX)
74    IF(ITYPE.NE.2.AND.ITYPE.NE.4) GO TO 84
      IF(XMN .GE. 0.) XMIN=IFIX(XMN)
      IF(XMN .LT. 0.) XMIN=IFIX(XMN)-1
      IF(XMX .GE. 0.) XMAX=IFIX(XMX)+1
      IF(XMX .LT. 0.) XMAX=IFIX(XMX)
84    A=YMAX-YMIN
      B=XMAX-XMIN
      IF(A .NE. 0.) GO TO 85
      YRANGE=.5
      YMIN=YMIN-1
      YMAX=YMAX+1
      GO TO 86
85    YRANGE=1./(YMAX-YMIN)
86    IF(B .NE. 0.) GO TO 87
      XRANGE=.5
      XMIN=XMIN-1
      XMAX=XMAX+1
      GO TO 88
87    XRANGE=1./(XMAX-XMIN)
C   COMPUTE INFO FOR AXES
88    XAXIS=(XMAX-XMIN)/FLOAT(ILEN)
      ILEN1=ILEN+1
      XLABEL(1)=XMIN
      DO 89 I=1,ILEN
89    XLABEL(I+1)=XLABEL(I)+XAXIS
      IF(ITYPE.NE.4.AND.ITYPE.NE.2.AND.ITYPE.NE.7.AND.ITYPE.NE.5)
     * GO TO 90
      DO 91 I=1,ILEN1
91    XLABEL(I)=10.**XLABEL(I)
90    XTOP=XRANGE*PLEN
      YLIM=6.*PHITE
      YLIM=IFIX(YLIM)
      YTOP=YLIM*YRANGE
      IYLIM=IFIX(YLIM)+1
      YAXIS=(YMAX-YMIN)/YLIM
      DO 92 I=1,IFIN,2
      WORK(I)=(WORK(I)-XMIN)*XTOP + 1.5
      WORK(I+1)=(WORK(I+1)-YMIN)*YTOP + 1.5
92    CONTINUE
C
C   START PLOTTING
C
      III=MIN0(6,MAX0(1,IFIX(ALOG10(YMAX)+2)))
      IF(ITYPE.EQ.3.OR.ITYPE.EQ.4.OR.ITYPE.EQ.6.OR.ITYPE.EQ.7)
     * III=MIN0(6,MAX0(1,IFIX(YMAX+2)))
      FORMAT(3)=CODES(III)
      FORMAT(6)=PLC(ILEN)
C   GET AXIS ANNOTATION INFO
      ILEN10=ILEN*10
      ICH=0
98    ICH=ICH+1
      IF(ITITX(ICH) .NE. 0) GO TO 98
      IXT=ICH-1
      ICH=0
97    ICH=ICH+1
      IF(ITITY(ICH) .NE. 0) GO TO 97
      IYT=ICH-1
      IF(IXT .EQ. 0) GO TO 95
      IXM=MIN0(IXT,ILEN10)
      ISP=(ILEN10-IXM+2)/2
      WRITE(IWRITE,803) (IBLK,I=1,ISP),(ITITX(I),I=1,IXM)
803   FORMAT(T9,122A1)
95    IYM=MIN0(IYT,IYLIM)
      ISY=(IYLIM-IYM)/2
C
C   PLOT X AXIS
C
      WRITE(IWRITE,801) (XLABEL(KKK),KKK=2,ILEN1,2)
801   FORMAT(T4,5(11X,1PE9.2))
      WRITE(IWRITE,800) (XLABEL(KKK),KKK=1,ILEN1,2)
800   FORMAT(T5,5(1PE9.2,11X),1PE9.2)
      DO 93 I=1,ILEN10
      LINE(I)=IMNS
      IF(MOD(I,10) .EQ. 0) LINE(I)=IPLS
93    CONTINUE
      ILEN11=ILEN10+1
      LINE(ILEN11)=IOH
      WRITE(IWRITE,802) (LINE(I),I=1,ILEN11)
802   FORMAT(T10,'O+',112A1)
C
C   DO ONE Y LINE AT A TIME
C
      ILEN12=ILEN11+1
      DO 100 III=1,IYLIM
      II=IYLIM-III+1
      DO 101 LL=1,ILEN11
101   LINE(LL)=IBLK
      M=II-1
      ITEST=MOD(M,6)
C   SCAN WORK ARRAY FOR VALUES OF Y EQUAL TO THIS Y VALUE
      DO 102 JJ=1,IFIN,2
      IY=WORK(JJ+1)
      IF(IY .NE. II) GO TO 102
      IX=WORK(JJ)
      JJ2=(JJ+1)/2
      DO 103 IP=2,IPOS
      IPS=IP-1
      IF(IND(IPS).LE.JJ2.AND.INUM(IPS).GE.JJ2) GO TO 104
103   CONTINUE
104   ICH=IC(IPS)
      IF(LINE(IX) .EQ. IBLK) GO TO 105
      IF(LINE(IX) .EQ. ICH) GO TO 102
      ICH=MIX
105   LINE(IX)=ICH
102   CONTINUE
      LINE(ILEN12)=IMNS
      ICH=IBLK
      IF(IYT.EQ.0 .OR. III.LE.ISY .OR. III.GT.IYT+ISY) GO TO 99
      ICH=ITITY(III-ISY)
99    IF(ITEST.EQ.0 .OR. III.EQ.1) GO TO 111
C   PRINT A LINE
      WRITE(IWRITE,900) ICH,(LINE(KKK),KKK=1,ILEN12)
900   FORMAT(1X,A1,T10,'-',112A1)
      GO TO 100
C   PRINT A LABELLED LINE
111   YY=YMIN+M*YAXIS
      IF(ITYPE.EQ.3.OR.ITYPE.EQ.4.OR.ITYPE.EQ.6.OR.ITYPE.EQ.7)
     * YY=10.**YY
      DO 110 IJ=1,ILEN11,10
      IF(LINE(IJ) .EQ. IBLK) LINE(IJ)=IPLS
110   CONTINUE
      LINE(ILEN12)=IPLS
      WRITE(IWRITE,FORMAT) ICH,YY,(LINE(KKK),KKK=1,ILEN12)
100   CONTINUE
C
C   PRINT X AXIS
C
      DO 96 I=1,ILEN10
      LINE(I)=IMNS
      IF(MOD(I,10) .EQ. 0) LINE(I)=IPLS
96    CONTINUE
      LINE(ILEN11)=IOH
      WRITE(IWRITE,802) (LINE(I),I=1,ILEN11)
      WRITE(IWRITE,800) (XLABEL(KKK),KKK=1,ILEN1,2)
      WRITE(IWRITE,801) (XLABEL(KKK),KKK=2,ILEN1,2)
      IF(IXT.NE.0) WRITE(6,803) (IBLK,I=1,ISP),(ITITX(I),I=1,IXM)
C
C  CLEAN UP PARAMETERS
C
      ICALL=0
      IPOS=1
      KJ=1
      IFIN=-1
      RETURN
      END
      SUBROUTINE PPSETU(IX,X,Y,ICAR,XLAB,YLAB,IPLOT)
      REAL X(100),Y(100),XLAB(80),YLAB(80)
      COMMON/UAB/WORK(3656)
      XH=0.
      XL=0.
      YL=0.
      YH=0.
      IF(IPLOT.EQ.-2)GO TO 200
      IF(IPLOT.EQ.0)GO TO 200
      IF(IPLOT.EQ.-1)GO TO 100
      YA=0.
      YL=Y(1)
      YH=Y(1)
      DO 10 I=1,IX
      IF(Y(I).LT.YL)YL=Y(I)
      IF(Y(I).GT.YH)YH=Y(I)
10    YA=YA+Y(I)
      YA=YA/IX
      IF(0.90*YA.LT.YL)YL=0.9*YA
      IF(1.10*YA.GT.YH)YH=1.10*YA
100   CALL PPINIT(1,6,XL,XH,YL,YH)
200   CALL PPSTOR(IX,X,Y,ICAR,WORK)
      IF(IPLOT.LT.0)GO TO 400
      WRITE(6,600)
600   FORMAT('1  ')
      CALL PPLOT(6.0,4.0,XLAB,YLAB,WORK)
400   CONTINUE
      RETURN
      END
      SUBROUTINE PPSTOR(NP,X,Y,ICHR,WORK)
C
C   STORAGE ROUTINE FOR PRINTER PLOTTING PACKAGE
C
C
C      WHERE:
C
C        NP - NUMBER OF POINTS TO STORE (FOR THIS PASS)
C         X - A VECTOR OF X VALUES
C         Y - A VECTOR OF Y VALUES
C      ICHR - THE PLOTTING CHARACTER
C      WORK - WORK VECTOR OF ATLEAST 2 TIMES THE TOTAL NUMBER
C                    OF POINTS STORED DURING THIS STAGE
C
C
C    PPSTOR MAY BE CALLED AS MANY AS 50 TIMES (TO OVERLAY
C    CURVES ON ONE ANOTHER) AFTER A CALL TO PPINIT.
C    A CALL TO PPLOT WILL PLOT THE DATA STORED BY THE PPSTOR
C    STAGE.
C
C
      DIMENSION X(1),Y(1),WORK(1)
      COMMON/PPBLK/ ITYPE,IWRITE,XMIN,XMAX,YMIN,YMAX,ICALL,IPOS,
     *              KJ,IFIN,IC(50),IND(50),INUM(50),FORCEX,FORCEY
      LOGICAL FORCEX,FORCEY
C
C   STORE THE VALUES FOR FUTURE PLOTTING
C
      N=NP
      IF(IPOS .EQ. 51) GO TO 60
      J=IFIN
      ISTART=ICALL
      DO 4 I=1,N
      IF(FORCEX .AND. (X(I).GT.XMAX .OR. X(I).LT.XMIN)) GO TO 4
      IF(FORCEY .AND. (Y(I).GT.YMAX .OR. Y(I).LT.YMIN)) GO TO 4
      J=J+2
      WORK(J)=X(I)
      WORK(J+1)=Y(I)
4     CONTINUE
      IF(J .EQ. IFIN) RETURN
      IFIN=J
      N=(J+1)/2-ICALL
C   CONVERT X AND Y VALUES IF NECESSARY
      GO TO(40,10,20,30,10,20,30),ITYPE
C
10    DO 11 IN=1,N
      I=2*(IN+ISTART)-1
      IF(WORK(I) .LE. 0.) GO TO 50
11    WORK(I)=ALOG10(WORK(I))
      GO TO 40
C
20    DO 21 IN=1,N
      I=2*(IN+ISTART)
      IF(WORK(I) .LE. 0.) GO TO 50
21    WORK(I)=ALOG10(WORK(I))
      GO TO 40
C
30    DO 31 IN=1,N
      I=2*(IN+ISTART)-1
      IF(WORK(I).LE.0. .OR. WORK(I+1).LE.0.) GO TO 50
      WORK(I)=ALOG10(WORK(I))
31    WORK(I+1)=ALOG10(WORK(I+1))
C
C   RECORD OTHER INFO
C
40    IC(IPOS)=ICHR
      IND(IPOS)=ISTART+1
      INUM(IPOS)=N+ISTART
      ICALL=ICALL+N
      IPOS=IPOS+1
      RETURN
C
C  ERROR MEAASGES
C
50    WRITE(IWRITE,700)
700   FORMAT(/' ***ERROR*** NEGATIVE OR ZERO VALUE ENTERED FOR '
     * ,'LOG PLOT.')
      GO TO 9
60    WRITE(IWRITE,701)
701   FORMAT(/' ***ERROR*** MAX PLOTS EXCEEDED - NO PLOT PRODUCED')
9     ICALL=0
      IPOS=1
      KJ=1
      IFIN=-1
      RETURN
      END
